General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



NASA CONTRACTOR REPORT 


166565 



Solution of Three-Dimensional Time-Dependent Viscous Flows 
Part 3: Application to Turbulent and Unsteady Flows 


(NASA-cs-ueses-pt-i) solution of 

3 -DIMENSIONAL TIME-DEPENDENT VISCOUS FLOWS- 

PAfiT 3: APPLICATION TO TUEBULENT AND 

UNSTEADY FLOWS (Scientilic Research 

Associates, Inc.) 93 p HC A05/MF A01 G3/34 


Bernard C. Weinberg 
Henry McDonald 



N84-27003 

Unclas 

13670 




NASA CONTRACTOR REPORT 


166565 


Solution of Three-Dimensional Time-Dependent Viscous Flows 
Part 3: Application to Turbulent and Unsteady Flows 


Bernard C. Weinberg 
Henry McDonald 

Scientific Research Associates, Inc. 
Glastonbury, Connecticut 


Prepared for 

Ames Research Center 

under Contract NAS2-10016 


IVJASA 

National Aeronautics and 
Space Administration 

Ames Research Center 

Moffett Field, California 94035 


TABLE OF CONTENTS 


Page 

SUMMARY 1 

INTRODUCTION . . . 2 

LIST OF SYMBOLS 4 

ANALYSIS 7 

Background 7 

Coordinate System 11 

Governing Equations 13 

Continuity Equation 14 

Momentum Equations 15 

Diffusion Terms ...... 17 

Convective Terms ..... 18 

Energy Equation 19 

Energy of State 20 

Linearizations 20 

The Turbulence Model ..... 21 

Spatial Difference Approximations ..... 23 

OR Operator Notation 23 

Application of Coupled Nonlinear Parabolic Equations .... 27 

Linearized Block Implicit Scheme . 28 

Implementation of the LBI Scheme Employing the OR 

Operator Technique 34 

Boundary Conditions and Initial Conditions 36 

The Computer Code 37 

Numerical Results 40 

Two-Dimensional Steady Turbulent Flow 40 

Weighardt Flat P'Late 40 

Newman Airfoil Case 44 

Two-Dimensional Unsteady Flow . 46 

CONCLUSIONS 51 


ACKNOWLEDGEMENT 


52 


TABLE OF CONTENTS (Cont.) 


Page 

APPENDICES 53 

Appendix A - Linearization Technique . . , 53 

Appendix B - Geometric Properties of Airfoils , , , . . 58 

Appendix C - Gens^Talized Operator Compact Implicit Schemes ... 70 J 

REFERENCES 74 

TABLES 

FIGURES ‘ 83 


M 

Ir 


SUMMARY 


The flow over a helicopter rotor is an Important example of three-dimensional 
time-dependent viscous flow. The boundary layers that develop on the rotor blades 
play a significant role in that they set loss levels and control retreating blade 
stall. Consequently, there is considerable interest in developing a numerical scheme 
for solving the time-dependent, three-dimensional compressible viscous flow equations 
in order to predict such flow fields, and which can be used as an aid in the design of 
helicopter rotors. 

In order to further investigate the numerical procedure, we have exercised the 
computer code that was developed under a previous phase of the current research 
program to solve an approximate form of the three-dimensional unsteady Navier-Stokes 
equations employing a Linearized Block Implicit technique in conjunction with a QR 
operator scheme. Results of calculations are presented for several two-dimensional 
boundary layer flows including steady turbuler*- and unsteady laminar cases. 

A comparison of fourth order and second order solutions indicate that Increased 
accuracy can be obtained without any significant increase in cost (run time). The 
results of the computations also indicate that the computer code can be applied to 
more complex flows such as those encountered on rotating airfoils. Finally, th 
geometry of a s}munetric NACA four digit airfoil is considered and the appropriate 
geometrical properties are computed. 
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INTRODUCTION 


The behavior of boundary layers on wings and bodies has long been of interest 
to aerodynamicists. In both steady and unsteady flows, the boundary layers are known 
to govern a major portion of the losses and to significantly influence the vehicle 
lift and moment coefficients. When the flow is steady, boundary layer prediction 
schemes based on numerical solution to the governing partial differential equations 
of motion have reached a high level of sophistication and predictive accuracy, even 
in three space dimensions. In unsteady flows, such as are commonly encountered in 
rotary winged aircraft, 'jome progress has been made in two space dimensions but 
little to date has appeared on unsteady three-dimensional boundary layers. 

Two particular problems arise with time-dependent three-dimensional boundary 
layers relative to the steady case. The first of these is the rather obvious one of 
time integration with its added requirements of transient accuracy coupled with an 
increase in the computational labor. The second of these is the so-called negative 
cross flow problem, which to some extent has troubled the steady boundary layer 
prediction schemes. Kendall, et al /’Ref. 1) discuss the negative cross flow problem 
for steady three-dimensional boundary layers in a very illuminating fashion. This 
particular problem arises when the spanwise component of velocity changes sign and 
will be discussed in detail subsequently. Because of the interest by external aerb- 
dynamicists in swept wing boundary layers where the negative cross flow problem 
(in this case flow from tip to root) is not usually encountered, the negative cross 
flow problem has not received a great deal of attention to date. However, in transient 
flows, particularly those encountered on rotor blades in forward flight, negative cross 
flows are frequently encountered. For instance, the advancing rotor blade has cross 
flows of one sign during the first ninety degrees of rotation and these can change sign 
over part of the blade during the second ninety degrees. 

Thus to be of practical value, time-dependent three-dimensional boundary layer 
prediction schemes require high computational efficiency and transient accuracy 
coupled to the ability to treat arbitrary cross flow profiles. 

In this report we describe the development of a computer code for the efficient 
solution of three-dimensional time-dependent viscous flows on fixed and rotary 
aircraft. The Linearized BJock Implicit (LBI) technique of Briley and McDonald (Ref. 2) 
in coordination with a tridiagonal QR operator scheme (Ref. 3) is employed to solve 
the reduced turbulent Navier-Stokes equations which are derived for nonorthogonal 


coordinates in generalized tensor form. The rationale for the choice of this 
approach is discussed in detail in Ref. 3 and 4. 

The basic assumptions made in the derivation of these equations are that the 
pressure does not vary normal to the shear layer and that in the energy equation 
the square of the normal velocity is neglected with respect to the other velocity 
components (To » constant) . The latter assumption is included only for computations 
simplification purposes and is not essential in the analysis. For turbulent flows, 
a two-layer mixing length model is employed and its formulation in generalized 
tensor notation is given. A novel method is employed for solving the continuity 
equation in conjunction with the reduced Navier-Stokes equations. The continuity 
equation is split by employing the Douglas-Gunn procedure to obtain a consistent 
approximation to the full equation which is then solved as an integral. Results of 
several boundary layer calculations are presented and comparisons with experimental 
data and other reported computations are made. 
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AI^ALYSIS 

Background 


Three-dimensional boundary layers occur on the wings and fuselages of both 
fixed and rotary wing aircraft, In both types of vehicles, the boundary layers are 
important in setting loss levels and determining useful operating ranges. As is 
well known, boundary layers are sensitive to pressure gradients. In time-dependent 
flow the temporal acceleration terms appear in the momentum equation in a form very 
similar to the conventional imposed pressure gradient and so for qualitative evalua- 
tion purposes can be regarded as "pseudo" or "auxiliary" pressure gradients, Viewed 
in this manner, the temporal acceleration terms can be seen to influence quantities 
of practical importance such as skin friction, displacement thickness and the onset 
of separation. At the range of frequencies typically encountered in rotary wing 
aircraft aerodynamic problems, it is clear, for instance, from the extensive review 
of McCroskey (Ref, 5), that very significant transient boundary layer effects can 
be observed. 

In examining the flow problems of practical Interest such as loss levels or the 
onset of separation it is evident that all three space dimensions must be considered. 
In conventional aircraft, the sweep effect is of interest and inherently three- 
dimensional. In rotary wing aircraft, in forward flight clearly very substantial 
transient changes occur in what might be termed the local sweep angle. However, 
generally speaking, the boundary layers remain thin unless catastrophic flow 
separation occurs or the flow at the wing or rotor tip is considered. 

Conventional boundary layer integration schemes have developed by forward 
marching the streamwise velocity u in the streamwise x direction and simultaneously 
marching out along the span in the 2 positive direction. In general, the spanwise 
marching scheme does not normally encounter negative w, i.e., spanwise inflow. This 
is very fortunate because it is difficult, indeed it could be argued impossible, to 
structure a physically satisfactory unconditionally stable noniterative scheme which 
permits forward marching in the spanwise direction with a negative w cross flow. 

At least Intuitively the problem of negative cross flow implies information being 
transferred upstream against the spanwise marching direction. Conventional stability 
analyses confirm the inability to forward march into regions of significant negative w 
From experience with attempts to march the two-dimensional boundary layer equations 
into a region of separated flow and its obvious relationship to the negative cross 


flow problem, it in not eurprising that spanwiso inarching into a negative croee 
flow regi n is not aceompliBhed without special treatment. Recently, conventional 
boundary layer developers have been turning to performing an implicit spanwise 
construction to remove the restriction of only positive cross flows (Kendall, et al, 

Rcf.l), tin and Rubin (Ref. 6) in their predictor-corrector boundary region 
solutions for flow over a y'iwed cone at moderate incidence have also shown that 
retaining diffusion in the spanwise direction not only eliminates the problems 
associated with negative cross flows, but improves upon the solutions obtained by 
standard three-dimensional boundary layer techniques. Furthermore, boundary 
conditions applied at the tip con influence the flow inboard, if required by the 
physics of th^^ flow. 

For these reasons, the implicit spanwise construction has been a feature of the 
three-dimensional duct flow analysis of Briley (Ref. 7) and McDonald and Briley 
(Ref. 8). As a consequence of these observations and the need to remove the negative 
cross flow restriction, a spanwise implicit formulation seems mandatory for 
rotary wing applications and ot least desirable for fix'^d wing applications, 
expeclally as it can be had for a very modest increase in code computational labor. 

Based on the experience in Refs. 7 and 8, the spanwise implicit sweep would only 
result in a moderate increase in computational effort relative to the explicit 
spanwise marching approach. The extension of the conventional three-diraesional 
boundary layer equations to allow spanwise diffusion is easily accomplished, and in 
view of the improved physical representation which thus follows, it is recommended 
and has been implemented in this effort. 

As a matter of course, it has been assumed that normal to the wall an implicit 
formulation would be structured. In recent years for boundary layer type problems 
there has been little dispute as to the efficiency gains to be had from an implicit 
formulation normal to the wall (Ref. 9). However, in the strearnwise direction for 
steady 2-D flow, the equations are normally forward marched and the implicit 
stability obtained entirely from being implicit in the normal to the wall direction. 

In time-dependent flows a similar structure is to be had so that at each time level 

one strearni'/ise (explicit) forward marching sweep could be made with two implicit 

sweeps in the spanwise and normal directions to give the desired unconditional » 

stability. As mentioned earlier, the explicit sweep would probably require less 

computational effort by about 20% than an implicit strearnwise sweep and of course 

less storage. However, since the solution is being time marched the opportunity to 

take a strearnwise implicit sweep at roughly the same cost as the explicit sweep 
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docs arise. If one does perform a streainwlse explicit sweep, then the linearisation 
of nonlinear terms is performed about the known spatial marching level . If an 
implicit irtreamwlse structure is adopted, then full time linearization can be 
utilized, Ttiat is the linearization of the nonlinear terms is performed about the 
known time level . As is pointed out in Ref. 8, it is easier to obtain a consistent 
high order accurate spatial- temporal linearization by marching in time than in space 
(in time the marching derivatives have the form u^^ whereas in space they are non- 
linear and have the form u^u^). Further, by structuring implicitly in the space 
marching direction, small regions of axial reverse flow would be permitted. As a 
result of the combined benefits of consistent high order linearizations and the 
inclusion of small separated zones, a strearawise implicit structure is advocated 
and has been employed in this effort, 

Since a full 3-R spatial integration is carried out at each time step of a 
transient calculation, spatial accuracy plays a very important role in the overall 
efficiency of the numerical method. By choosing a spatial scheme of sufficient 
accuracy, one can expect to reduce the total number of grid points to an acceptable 
level. However, for the chosen scheme one must account not only for the spatial 
discretization errors but those introduced as a result of the linearization procedure. 

In order to get most out of a given spatial difference formula, the errors 
incurred in repr?' 'XJ-lag nonlinear terms by linear combinations of terms should be 
less than ar equal to the spatial discretization errors. If the linearization intro- 
duces a greater error than the spatial differencing, then either a coarser spatial 
mesh could be used, or Iteration, or some form of linearization improvement is called 
for. Iteration across a time step is not recommended since this only reduces the 
linearization error and computationally costs as much as a complete time step. 

Cutting back the time step would be prefercable to Iterating to preserve the lineariza- 
tion error at some acceptable level, since cutting back on the time step would Improve 
both the transient error and the linearization error. This point is clearly demonstrated 
in Ref, 3. To obtain a linearization, which introduces errors of at most the same 
as the spatial difference formulae, a Taylor series expansion about the known time 
level can be performed. This process clearly demands a formal block, i.e., coupled, 
treatment of the system of equations. For instance, in the streamwise momentum 
equation a typical term is linearized: 

(uw) «u w +uw -uw+0( At*) 
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n+1 

ond clearly one cannot lag w at the old time level n without introducing a first- 
order time error in order to get an uncoupled system, i,e», w not appearing in 
the stroarawise momentum equationp Thus, formal linearization and consideration of the 
resulting errors indicate the coupled system ought to be treated from the accuracy 
point of view. This is further reinforced when it Is realized that block, i.e., 
coupled, systems are not computationally expensive (in a relative sense). 

Additionally, a second type of approximation arises unconnected with lineariza- 
tion but arising from basic coupling terms in the original equations and if indeed 
some term:; in an equation are time lagged in order to uncouple the equation system 
and these terms are of equal importance to the terms retained, then again an Iterative 
updating is called for in order to achieve stability, accuracy and consistency. 

This could be termed ad hoe equation uncoupling, Blottner (Ref. 9) has shown that 
many iterations around the ^ hoc uncoupled set (> 10) are sometimes required in 
order to achieve an overall solution accuracy commensurate with the local difference 
molacule accuracy. 

The linearization technique discussed above is described in Ref. 8, together 
with its application to block coupled splitting schemes. Schemes of this general 
type are here referred to as "split linearized block implicit" or split LBI schemes, 
and are reviewed in detail by Briley and McDonald (Ref. 2). 

As a general observation, care is required to obtain acceptable transient 
accuracy for long time integration with conventional finite difference schemes. 

A Crank-Nlcolson centered time-implicit scheme for instance, although second order in 
time, shows quite a dispersion problem (relative to other schemes) on the simple 
pure convection problem. However, the problem of transient accuracy is significantly 
reduced in the typical boundary layer problem since the time dependency is con- 
tinuously input through boundary conditions and the concern is with relatively 
"short" time integrations. The computational problem is more of what the phase lag 
of the wall shear Is, relative to the prescribed free stream disturbance, than 
concern over the convection velocity of a wave in a shear after a long propagation 
time. The interest is in forced oscillations with a minimum scale of the boundary 
layer thickness over a few cycles of the motion, just enough to obtain repetition 
cyclically. It is, therefore, expected that a significant dispersion problem will 
not arise with a conventional implicit scheme. 

The governing equations that are considered here are the Navier-Stokes ei-uations, 
continuity, energy and the equation of state which are written in generalized tensor 
form for a body oriented coordinate system (boundary layer coordinates) . In 


10 


accordance with the boundary layer assumptions, the normal momentum equation is 
eliminated and the pressure is specified throughout the viscous layer in its stead. 

For the energy equation, constant stagnation temperature is assumed. This 
assumption is a good appioximation for the flow fields considered, and is thus 
included here only for purposes of simplification. In the analysis that follows, 
the full energy equation could equally well have been used. Employing the equation 
of state which relates the pressure p to the velocity components u and w by an 
algebraic equation, the problem can be reduced to one involving only the three 
velocity components, u, w and v and three equations, the streamwise and spanwise 
momentum equations and the continuity equation. Hence, we consider a block-three 
system rather than a block-four system which leads to a significant reduction in 
computer time. If the full energy equation were to be considered, a block-four 
system would result due to the inclusion of the temperature as an additional unkno\\rn. 

Coordinate System 

Since the goal of this effort, is to solve for the flow over airfoils an under- 
standing of the type of geometries to be considered is essential to guide the 
choice of the coordinate system and the structure of the computer code. The coordi- 
nate system is not only dependent upon the geometry of the airfoil, but also upon the 
approximations that are made to the governing Navler-Stokes equations, ‘As in boundary 
layer theory, we also assume that in the approximate form of the Navier-Stokes 
equations the pressure is constant normal to the shear layer. Inherent in the 
assumptions is that the shear layer is thin. As pointed out by Howarth (Ref. 10) 
the boundary layer assumptions lead to the conditions that one coordinate directio'n 
must be normal to the body surface (this being straight) while the other coordinate 
directions must lie on the body surface. These conditions uncouple the metric data 
on the surface from that in the normal direction. Hence, the metric data for the 
surface coordinates are functions of the surface coordinates alone, while the metric 
data for the normal coordinate direction are functions of that coordinate alone. 

The choice of the surface coordinates is rather arbitrary and is based on considera- 
tions such as the ease of construction or the grid distribution on the wing surface. 

In the numerical solution of the flow over an airfoil there are many advantages 
to be gained by the judicious choice of coordinates. The most obvious advantage is 
that the physical boundaries of a flow region can be represented by coordinate 
surfaces, '^his removes the need for fractional cells in general; hence, the compli- 
cations and loss of accuracy associated with a boundary interpolation are removed. 
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Another advantage is that a uniform numerical mf'f'od can be used. The solution can 
then be performed with a fixed number of cells in any given direction and with a 
uniform mesh spacing. For a rectangular plan form wing a Cartesian coordinate 
system would be adequate. However, for more general plan forms such an c swept wing 
a nonorthogonal grid which conforms with the boundaries is preferred since it covers 
the entire airfoil while a Cartesian grid would not. 

Another consideration is the selection of a coordinate grid distribution; 
the major objective being the resolution of large solution gradients. The approach 
taken here is to construct coordinate transformations that contain distributions for 
physical nesh points. In this context, ths uniforro mesh of computational space is 
simply mapped into a suitably distributed mesh in physical space. When the trans- 
formation contains the mesh point distribution, there is no need to construct the 
apparatus for the discrete approximation of derivatives on a non'iniforra mesh. This 
results in a savings in both computer logic and storage. 

Hence, in this work, a coordinate system is chosen that conforms with the 
boundaries of the physical domain, i.e., the wing surface which in general will be 
nonorthogonal. In addition, provisions are made for analytical grid transformations 
(Ref. 11) in each coordinate direction, in order to suitably distribute grid points 
in regions of large gradients. 

In view of the type of geometries to be considered and the assumptions made to 
obtain the approximate form of the Navier-Stokes equations, a specialized nonorthogonal 
coordinate system is advocated where the metric tensor which has four independent 
components is given by 


^11 

9 12 

0 

^12 

9 22 

0 

0 

0 

933 


The subscripts 1 and 2 refer to the directions on the surface of the body while 
subscript 3 refers to the direction normal to the body. Since the metric data in 
the coordinate directions on the airfoil surface are not functions of the normal 
direction, the metric data in a 1 - 2 surface above the body are evaluated on the 
body surface (Ref. 10). Furthermore, due to the use of nonorthogonal coordinates 
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it becomes advantageous to derive the equations In general nonorthogonal coordinates 
employing generalized tensors. Details of the generalized tensor notation can be 
found in Refs, 12, 13, and 4, 

An important feature of the analysis to follow is that the governing equations 
which are derived, under the prescribed assumptions, are invariant for any coordinate 
system or any grid transformation (although, of course, the physical approximations 
are coordinate dependent). The grid transformations are absorbed into the geometrical 
coefficients, leaving the equations unaltered in form. This feature allows for the 
construction of the geometric data to be contained in one subroutine with the defini- 
tion of the metric data and their derivatives as input. 

Other forms of the governing equations for nonorthogonal coordinates which do 
not rely on generalized tensor notation were considered, i.e,, such as the steady 
three-dimensional boundary layer equations as given by Cebeci (Ref. 14 ). Since in 
that set of equations the only viscous terras retained are these in the normal direction, 
any nonorthogonal effects are Introduced solely through the convective terras. This 
is a direct consequence of the assumptions which decouple the normal coordinate 
from the surface coordinates. In those equations, the geometric terms that explicitly 
appear are the metrics, h^^ and h^, the geodesic curvatures kg^ and kg^, the angle 
between the surface coordinate lines and an additional nonorthogonal curvature terra. 

In contrast to Cebeci *s equations, the set of equations considered here, the 
approximate form of the Navier-Stokes equations, allow for all the viscous terms to 
be retained. A description of these terms without generalized tensor notation would 
be curabersome. Hence, the generalized tensor approach was chosen. 

In Appendix B the geometric properties of a surface in three dimensions are 
discussed and where appropriate, the generalized tensor equivalents are given. 

In addition, a symmetric (uncarabered) NACA four-digit airfoil is considered and the 
pertinent geometric coefficients are presented. 

Governing Equations 


In view of the ultimate goal of this program, to solve an approximate form 
of the unsteady three-dimensional Navier-Stokes equations on airfoil shapes, the 
governing equations are derived in general nonorthogonal coordinates and are 
given in generalized tensor notation. 

In the following derivation the governing equations are nondimensionalized as 

follows, x^ with respect to the characteristic length L, the velocity with respect to 

2 2 

U^, density, pressure and temperature with respect to p^, and /c^ respectively 

and time with respect to L/U^. Viscosity is nondimensionalized with respect to 


13 



Continuity Equation 

Consider the continuity equation wlrtten in vector form so that it is indepen- 
dent of coordinate system, i.e., 


at 


+ V • ^ q =0 


( 1 ) 


where p is the density and q is the velocity vector , which can be written 
in a covariant basis as 


( 2 ) 


where u is the 1-th contravariant velocity component and e. is the covariant basis 
vector in the x direction. The velocity vector could have been expressed in a number of 
different forms, each with its own attribute. Here for the moment the velocity vector 
is expressed in a covariant basis, for simplicity. Subsequently, the velocity vector 
will be transformed into its physical components for the numerical solution of the 
governing equations. The reason for this is that the contravariant basis exhibits 
variation in its components if the coordinates are such that the metric varies. 

For boundary layer flows, the physical velocity components are roughly aligned with 
the coordinates and exhibit no variation with the metric per se. As such, it is 
felt that the actual computations are better performed on the pnysical components. 

The divergence of a vector is given by 


■pq-- (pu'^) +pu' C = -y 


k. 1 k 

where pu is the covariant derivative, pu 

ih 


Ik ^ " ”k 

is the partial derivative in the x 


(3) 

k 


direction, J is the Jaccbian and is the Christoff el syitbol. 

In Equation 3 two forms of the divergence are presented, one involving the 
Chrlstoffel symbol or curvature term directly and the other the Jacobian. The 
former is perhaps more restrictive since it requires additional smoothness of the 
geometrical quantities. However, both forms are used depending on which is more 
convenient for the given application. For the momentum equations the form involving 
the Chrlstoffel symbols is employed in the evaluation of the explicit (lagged) dif- 
fusion terms while for the continuity equation the form involving the Jacobian is 
used. Thus, the continuity equation can be expressed as 


at 


(J^u ),j^=0 


(A) 
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Momentum Equations 

The momentum equations in vector form can be wrxtten as 


-p[ 




(q -Vjq ] 


where a is the stress tensor. 

In generalized tensor notation, Equation 5 becomes 


“‘‘“'Ik] «l = “-'Ik®! 


The stress tensor is defined as 


2 ^a1r"'+ -t- .IK 
3^AJS + < 


where y is the viscosity, p is the pressure, A is the velocity divergence and g 
is the strain tensor, and the Reynolds number. Re, is defined as P„U^b/y^, The 
strain tensor expressed in terms of velocity components is 

c"' = u' I g"”‘ + u"| g'"' (8) 

•m 'm 

where and are the components of the metric tensor. Employing the fact that 
and substituting the definition of the strain into the stress tensor, 

we obtain 

(P T ^ W H 9 + ^ “ l„9 (9) 


Substituting Equation 9 into the momentum equation and employing the relationship 

g*'^| = 0 (10) 

'k 

we obtain for the i-th momentum equation in the e^ direction 


r k il 1 

9|kJ 

= -g"'(p + 

2 

3 

-^").k 


(11) 
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In Ref, 3 it was pointed out that the QR operator scheme requires that derivatives 
in any direction operate on only one variable. In the momentum equation this require- 
ment prevents the Implicit treatment of certain diffusion terms that arise due to the 
curvature of the body, Although these terms are often treated explicitly anyway the 
use of standard finite difference techniques instead of the QR operators would give 

one the opportunity to treat these terms implicitly, if so desired, This and the 

use of the quasi-llnear form of the governing equations are the only major limita- 
tions that arise in the QR operator treatment of the approximate form of the 

Navler-Stokes equations. In the usual boundary layer approximations, these explicitly 

- 1/2 

treated terras would not appear in the equations since they are of order 0 (Re ' ) or 
smaller, and should, therefore, be of little consequence. In principle, the quasi- 
linear and, for instance, the full conservative for-' of the differential equations, 

are equivalent. In discrete form, various formulations of the governing equations 

exhibit different properties (Ref. 15). In the present problem, no distinct dis- 
advantage appears to arise from the required use of a quasi-linear form of the 
governing equations. 

The requirement that derivatives in any direction operate only on one variable 
would be more restrictive in the treatment of the 'pressure gradient term in the full 
Navier-Stokes equations, The linearization of this term introduces derivatives of 
all the velocity components in a given direction. According to the limitations of 
the QR operator scheme described above, some of these terms must be treated explicitly. 
Since an explicit treatment of these terms could reduce the stability bound of the 
calculation scheme, an alternate procedure should be considered. This would involve 
the addition of an auxiliary equation relating the pressure gradient term to the 
derivatives of the velocity components and would increase the block size of the 
system. An assessment of the efficiency of such a procedure has not been carried 
out and further work in this area would appear to be warranted. 

In the discussion that follows, the governing equations are first split into 
an explicit part and an implicit part in accordance with the QR operator requirements. 
Thereafter, the resulting equations are cast into "standard form", so that the 
equations can be appropriately linearized and treated with the LBI technique. 

Since mixed partial derivatives are commonly treated explicitly in orthogonal 
coordinate systems, a similar procedure will be used for generalized nonorthogonal 
coordinates where the explicit treatment is extended to include mixe d second covariant 
derivatives. All other second covariant derivatives are retained as implicit. 

Although such a procedure would automatically treat more terms explicitly than one does 
for orthogonal coordinates, it simplifies the bookkeeping requirements in the 
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construction of the computer code, and is thus adopted here. Conversely, retaining 
the covariant derivative as a unit (whether implicit or explicit) could prevent 
instabilities that may arise due to time splitting. This can occur when two 
portions of one term should cancel identically, but do not due to their being 
split between two sweeps. 

Diffusion Terms 

Consider the term 

( 12 ) 

If j = k the term is retained as implicit, and if j k then it is treated explicitly. 
We will consider the case j / k first. Upon expanding the explicit part of the 
diffusion term, it becomes 



Note that the first term on the right-hand side of the equation is in conservation form. 
Although the implicit equations are treated in quasi-linear form, for the purpose of 
evaluating the explicit terms the most convenient representation is used. Tfje 
implicit terras, with j = k become (note there is no sum on j) 




where 


and 



r 


m 


jn * mj 


- r 


nm 

mn * jj 


+ r' 


jn,j 




U 


m pi 

•j *mj 




(15) 


( 16 ) 


m=i+l, n=5-2i and there is no sura on m and n 

Since T^^ involves velocity compo'.ients and derivatives in directions other than the 
i-th direction, the term is also treated explicitly. 

Hence the total diffusion terras for the 1-th momentum equation is given in 
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quasi-Iinear form as 


§^fQfNAL PAGE fS 

Op Poor quality 



where 


and 


and repeated indicies do not 



0 0 “ 

I 0 

0 I > 

1 I " 

2 I 

I 2 - 


indicate summation. 


(18) 


(19) 


Convective Terms 


The convective term for the i-th equation can be written as 

'j ' L j rnj J 

j i I ^ f j i |. J m i " i 1 

pu u |j- I { />u u ..^•z py> u r^j+ &,ipu u r„j) 


which becomes upon expansion 

3 


( 20 ) 


( 21 ) 


where the last term is nonzero only when j = 3. The full momentum equation is 
obtained by substituting Equations (17) and (21) into Equation (11) and treating 
the pressure gradient and velocity divergence as explicit terras. Since the pressure 
is specified and impressed upon the viscous layer, its specification replaces the 
normal momentum equation. Thus, the streamwise and spanwise momentum equations are 
the only two retained. 


The momentum equation now beqomes 
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I { [c| [c| - 1 cJgVriJu'.j 


+ s|[|cyV.Krkl]u +[|ctg'Vsyu 


k kk J 1 J 

n 1 1 C > . . I 1 1 ' — 




i nl r k 


+ (l- S|[)g"[^u''l ] } 


inerEv Equation 


The energy equation employed here states that the stagnation temperature is 
constant throughout 

:k ( 23 ) 

2 

The generalized tensor notation q is given by 

q2 = u’yjgj. 

where u^ and u^ are the contravariant velocity components. Incorporating the assump- 
tions made concerning the coordinate system employed, i.e. 


g = g = 0 

13 23 


we obtain 


q 2 = (u')^ g,, + 2u'u^g,g + (u®) ^ g22 + (u®)^ g- 


Neglecting the term involving with respect to the other terms, and defining 

physical velocity components, i.e. 


Up = u'h, 


Wp = u^hg 
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we obtoin 


T„ = T + 


This is the form of the energy equation used. 

Equation of State 

The equation of state assumes a perfect gas and is given by 

_ r-> 

r 

Linearizations 


p= -y-^T 


(24) 


(25) 


The following analyses assume a set of linear partial differential equations. 
However, the continuity equation and the convective part of the momentum equation 
are nonlinear, containing terms that involve the product of density and velocity 
components. In order i:o overcome this difficulty we employ the linearization procedure 
(described in Ref. 8 and reviewed in Appendix A) to linearize the aforementioned terms 
by Taylor series expansion about the known time level solution. 

The density is first eliminated by employing the equations of state and energy, 
and thereafter the resulting terms are linearized, These terms are of the following 
form 






^12 /..2^n /..h" 




(26) 




h,h'2 


n +/3 


■ n 


,f>+^ 




where all velocity components are the contravariant ones, and 6 is always a velocity 
12 3 

component, (u , u , u ) while tjj can be either a velocity component or a derivative of 
a velocity component. In the case when only one velocity term appears, i.e. , pu^, 
we set ij/ = ijj =1. 
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m, 5 . 




Ic is important to note that in the preceeding equations the contravariant 
velocity components are used# However, as noted in Ref, lA, it appears advantageous 
to solve for the physical velocity components, Therefore, when the governing 
equations are subsequently cast into a form amenable to the application of the LBI 
scheme, they are transformed so that the physical velocity components appear. 

The Turbulence Model 

We treat the set of three-dimensional ensemble-averaged turbulent reduced 
Navler-Stokes equations. Ensemble averaging permits the appearance of low frequency 
(relative to the turbulence) time-dependent "mean" flow. It is, therefore, necessary 
to specify a turbulence model suitable for this problem. 

The approach taken in the present effort assumes an isotropic turbulent viscosity, 
py, relating the Reynolds’ stress tensor to mean flow gradients. 

Reynolds stress = 

Using Favre averaging (Ref, 16), the governing equations then are identical to the 
laminar equations with velocity and density being taken as mean variables and vis- 
cosity being taken as the sum of the molecular viscosity, y, and the turbulent 
viscosity, 

At this point we require additional closure assumptions for the Reynolds 
stresses, i.e., the evaluation of y^. There are a variety of approaches available, 
from the simpler mixing length models to the more complicated one and two-equation 
models. Since the intention here is to verify the code's performance in wall bounded 
cases, we chose the mixing length models which have worked well in the past for similar 
flow environments (Ref. 17 ). The extension to more complex models could be undertaken 
at a later time if warra.ited. At that time, the LBI procedure that is used for the 
solution of the momentum equation could be applied to the k and e equations. 

Employing the Prandtl mixing length concept, the turbulent viscosity is 
given as . 

A'-T ~ pi (28) 

where i is the mixing length, and e* is the mean strain, which for the general case 
is related to the dissipation function. 

e* = <29) 
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where 4 / is the dissipation function, In Cartesian tensor notation <J* is given as 


where 




p - ^ r iiiL + Ml 
*^1) - 2 I axj 7^, J 

In the case of two dimensions, does not reduce directly to I’g^l (where Xg is nornml 
to the wall) without the additional assumption that the other components arc small. 

In nonorthogonal coordinates, with generalized tensor notation employed, an 
analagous CApression can be obtained for the dissipation function, i.e, 


where 


e ^ U] . 


(J Im 1 , Jm I 

e = g u' + g u 

m m 


After some algebra, <l> reduces to 


e ■' e. 


, . I 

J I > Jm I 

e,. = g, , uM + g u 
ij !m m 

For computation purposes, 4> can be expressed somewhat simpler, i,e. 


+g^^u'| ]gj. 

^ ’m *m 


On expanding, and recalling that « 0 m 5 ^ 3, 4> can be expressed as a summation 


u'L “1i + g‘"u'|)+g”!:l:g, u'| u"| +g lig‘v| 

k j 'I' 'I k i r I m Ik 'm k I Ik 's '5 «m j 'm 


+ g„g^^u^ u® 

^ 33 ^ 33 
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As In tho Cartesian tensor formulation, 4- does not automatically reduce to the 
dominant shear for two-dimensional boundary layers, Hence, provisions were made In 
the computer code that on option retains only the dominant component of the strain. 


The mixing length formulation is based on the two-layer model 
inner layer 
outer layer 

where y is the normal distance from the wall, is the van Driest damping term 


(36) 


^ - t- expl-yVA’*') 


and the von Karman constant k « and the van Driest damping coefficient A ■ 26,0, 
the dimensionless distance y+ is defined as 


and the friction velocity u^ is taken as 


in as 


In the outer layer f is given by 


lcQ“ XS 


where 6 is the local boundary layer thickness defined as .995 u^ and \ .09. 

This form of the mixing length has a discontinuity in its first derivative 
at the matching point between and An alternate form which avoids the dis- 

continuity and varies smoothly between the two layers has been proposed by 
McDonald (Ref. 18 ), and is given by 

I s i^^janh [-J- ) 2) ^37 J 

^CD 

In the turbulent calculations to be presented we e ^loy Eq. (37). 


Spatial Difference Approximations 

QR Operator Notation 

In this section, Implicit tridiagonal finite difference approximations to the 
first and second derivatives and to the spatial differential operator are considered. 
The QR operator procedure for generating a variety of spatial discretizations is also 
introduced. As special cases, standard second-order finite differences, first- 
order upwind differences, fourth-order operator compact implicit (OGI), fourth- 
order generalized OCI and exponential type methods are obtained. Since all these 
schemes are of the same form (cf. below), a single subroutine which defines the 
difference weights is all that is required to Identify the method, while leaving 
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the basic structure of the program unaltered. Subsequently, the results of numerical 
experiments employing some of these schemes will be presented. The rationale for the 
use of the QR approach in the present problem is discussed in dctpxx in Ref, 3. 

The QR formulation allows for ADI methods and permits the treatment of systems of 
coupled equations, l.e,, LBI methods. Although variable mesh schemes can be employed 
within the QR framevzork, it is believed preterable to use analytic transformations to 
obtain a uniform computational mesh, hence attention is restricted to uniform meoh 
formulations. 

The general concepts and notation will be introduced for two-point boundary 
value problems and then the methodology will be extended to more gener<il linear and 
nonlinear parabolic partial differential equations in one dimension. The application 
of the QR operator method to multidimensional problems is discussed in the section 
pertaining to the LBI scheme. 

Consider the two-point boundary value problem 

U(u) »a(x)Uj^j^ + b(x)Ujj + c(x)u « f(x) (38) 

with boundary values u(0) and u(l) prescribed, Derivative boundary conditions, 

although not discussed here, can easily be incorporated into the framework of the 

QR operator notation. Let the domain be discretized so that x^ » (j-l)h, j =1, 2, 

. . . , J + 1, and U, ~ u(x.), F, ~ u (x, ), S. u (x.) and h = l/j is the mesh 

J J J X J J XX 2 

width. The numbering convention was chosen here to be compatible with FORTRAN coding. 

Without loss in generality for a(x) ^ 0, Eq. (38) can be divided by a(x) so that 
we may treat instead the following equation 

L(u) - + b( x)Ujj + c(x)u » f(x) (39) 


where ~ . 

b(x) - b(x) / q(x), c(x) • c(x) /a(x) ond f(x) ■ f(x)/a(x) 

Substituting the finite difference approximations to the first and second 
derivatives 




into Eq. (39) and rearranging, we obtain 


Uj.,-2Uj+Uj^l 




(40) 

(41) 


L(u) ~ Sj + bjFj 




(42) 
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or 


[l- +[h"c,-2]u) +[ 


I + 


Rc 


U 


j + l 


h^fi 


(43) 


where Rc^ « hb^ is the cell Reynolds number, 

Equation (43) can be generalized by introducing operator format, i.e., 

rj' U^., + rj Uj + r[ + qjfj + q| f (4^) 

where the superscripts (-) minus, (c) center, and (•{-) plus indicate the difference 
weight that multiplies the variable evaluated at the (j-1), (j) and (j+1) grid points 
respectively, and where the and qj's for grid point j are functions of h, 

^i“l’ ^j+1* ^j-1’ ^j+1’ Comparing Eqs. (43) and (44) we can identify 

the r^'s and qj's, viz., 


r" - I - RCj /2 


- I + Rcj/2 


qj » 0 


c 


Qj * 0 


(45) 


We now define the tridiagonal difference operators Q and R 


r[uj] 


rj Uj.| + r' Uj + r* u,„ 

°[fj] 


(46) 


Noting that L(u) = f and substituting Eq. (46) into Eq. (44), we obtain 

r[uj] .h’'o[»,]= h^o[L(u),] (47) 


Alternatively by employing the inverse operator Q ^ and expression for b(u)^ can be 
obtained , 


‘-(“>1 ' -p-O-'RUj 


(48) 
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For standard central finite differences Q * Q ^ « I, the identity matrix, 

(the spatial operator is given explicitly in terms of and so that 

nothing would appera to have been gained in obtaining Eq. (48). However, in general, 

-1 

for higher order methods whereas Q is tridiagonal Q is a full matrix. Hence, 

Eq. (48) provides us with a convenient expression for the spatial operator for a 
wider class of difference approximations. The formalism in Eq. (48) is also applicable 
for first and second derivatives appearing alone (cf. Ref, 19). It should be pointed 
out, however, that Eq. (48) is not the most general formulation since the compact 
implicit formulas cannot be combined to yield a single scalar equation relating the 
spatial operator to the function values (Ref, 19). 

In Refs. 3 and 20 a technique due to Berger, et al is described for constructing 
fourth order tridiagonal methods which pos.'^ess a monotonicity property as the cell 
Reynolds number is Increased, R^ A brief description of the method for 

deriving generalized OCI schemes is given in Appendix C and the resulting Q and R 
coefficients are presented in Table II, 

Another family of schemes that can be expressed in QR operator notation are the 
so-called exponential methods. The idea, originally due to Allen (Ref. 21) 
(independently derived by II 'in (Ref. 22) and McDonald (Ref. 23)), and employed by 
Dennis (Ref. 24) is to set the difference weights so that the numerical solvition is 
equated to the analytic solution for the locally frozen constant coefficient 
equation. The QR coefficients of this exponential scheme is given in Table III, 

This method is second order accurate for Rc = 0(1) and becomes first order accurate 
as Rc ->- <» where the scheme reverts to first order upwind differencing. 

Another exponential scheme which is uniformly second order accurate was 

developed by El-Mistikawy and Werle (Refs. 25 and 26). The "exponential box scheme" 

which is incorporated in their solution of the boundary layer equations with strong 

blowing, is based on a spatial operator of the form given in Eq. (39). Berger, et al 

(Ref. 27) derived the counterpart for an operator of the form given in Eq. (30), 

but with c = 0. The Q and R coefficients are presented in Table IV. Although this 

scheme reverts to second order upwind differences as Rc it does not possess a ' 

discrete maximum principle while the exponential scheme of Allen (Ref. 20) does. 

In Table V a centered finite difference scheme is presented which permits the 

addition of artificial dissipation to the spatial operator when i 1 that 

Irc! never exceeds Rc 
' ' max 
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Annllcatlon to Coupled Nonlinear Parabolic Equations 

Before considering the LBI technique, we discuss some of the limitations placed 
on the QR operator scheme in solving a system of nonlinear parabolic equations. 

Given a system of m nonlinear parabolic equations in m unknowns, 


m 

I 




, n+i n. 
At 


- Nj X|, XgjXj,!) 


I - 1,2,..., J+l 


n*4’6 

where is a quasilinear spatial operator, the QR formalism carries directly over 

provided that for any equation only one independent variable is operated upon by the 
differential operator, For example, 

I 

0(u,w,vT “l ■ “xx + b(u,v,w)u, + c(u,v,w) 

is allowed since x derivatives of u only appear, while 

I 

■a{L;w,v) + b(u,v,w)u^ + c(u,v,w) + d(u,v,w)w^ i 

is not allowed since x derivatives of both u and w appear. The approximate form of 
unsteady Navier-Stokes equations used here, when written in quasi-llnear form, falls 
within the class of allowable differential operators. Thus, for the problem being 
addressed in the present study, the OCI schemes are applicable. Note that within the 
splitting approach, non allowable terms in the OCI scheme such as dw above, may be 
split off and treated by a special implicit sweep. Provided care is taken and for 
instance the Douglas-Gunn formalism is adhered to, no particular problem arises 
other than the cost of an additional implicit sweep which is Incurred. 

Thus, multidimensional problems and/or more general equation forms can usually 
be accomodated by a splitting procedure, which reduces the differential operator to a 
sequence of one-dimensional problems which have the appropriate allowable form. 
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However, as with standard finite differences, to avoid the cost of additional 
implicit sweeps, special procedures must be applied to cross derivative terms, e.g., 
extrapolation or explicit treatment. 

Linearized Block Implicit Scheme 

Consider a system of nonlinear partial differential equations 

(49) 

where ^ is a vector of unknowns and 1* is a source term vector which is a function of 
x^, x^, x^ and t. Extension to source terms which are functions of ^ are discussed 
in Ref. (8). ^ is a three-dimensional nonlinear differential operator and the 
matrix A appearing in the momentum equations is equal to pi where p is the density 
and I the unity matrix. 

Equation (49) may be centered abut the n+B time level, i.e. t^*^^ = (n+B)At = 
nAt+BAt = t’^+BAt, and written 

( 50 ) 

where 0 - B - 1 is a parameter allowing one to center the time step, i.e., B = 0 
corresponds to a forward difference, B = 1/2 to Crank-Nicolson and B = 1 to a backward 
difference. 

til 

After linearizing Equation (50) by Taylor series expansion in time about the n 
time level by the procedure described in Ref. 8 to give a second-order linearization, 
we obtain 

-2."!% f'’*® ( 5 „ 

where^^is the linearized differential operator obtained from ^ by Taylor Series 
expansion in time. 

The difference between the nonlinear operator T> and the linear operator ^ is 
defined as ~ At the intermediate level n + B, is represented as 

Using these relationships and dropping the vector superbar for convenience a two- 
level hybrid Implicit-explicit scheme is obtained 
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a"(<3>"‘‘''-<3>")/AI = + m" 43 " + (53) 

n*4*6 

The vector t ^ represents all of the terras in the system of equations which are 
treated explicitly. More about this will be said later, but for the moment note that 
may be approximated to the requisite order of accuracy by some multilevel linear 
explicit relationship, or approximated by with a consequent order reduction in 
temporal accuracy. 

The I -era tor is now expressed as a sum of convenient, easily invertible sub- 
operators • ' * -^m‘ In the usual ADI framework these suboperators 

are associated with a specific coordinate direction. Further, it is supposed that 
these suboperators can be expressed in the QR notation introduced earlier. Writing 
and mV as a single source terra Equation (45) is written as 

To solve this system efficiently it is split into a sequence of easily invertible 
operations following a generalization of the procedure of Douglas and Gunn (Ref. 19) 
in its natural extension to systems of partial differential equations. The Douglas- 
Gunn splitting of Eq. (54) is written as the following three-step procedure 


a" [$•”$'’] /a, = -$"] 

+ [j” 
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which can be transformed to the alternative form 


[A"-At/3i3"][4>*“"-<J>"] = a"[<J>»*-<I>'’] 


( 56 ) 


If the intermediate levels are eliminated, the scheme can be written in the so-called 
factored form 

« (57) 

4>"+Ats"*^ 

At this point it becomes necessary to consider the structure of the operators 
-^1’ 2 It will be recalled from the one-dimensional scalar problem 

that use of the QR format greatly facilitated the introduction of a wide variety of 
spatial difference fernmlae. It follows that in the extention to multidimensions 
undertaken here, the use of the QR formulation results in the appearance of the 
inverse operator with the sub-blocks of the ^ 3 oP®^ 3 bors. In 

order to implement the scheme the inverse operator Q ^ must be cleared. Accordingly, 
the scalar operator Q is generalized to the vector operator with (diagonal) sub- 
blocks In this generalization j = 1, 2 apply to the momentum equations and j = 3 

applies to the continuity equation. The i subscript is associated with the coordinate 
directions of the ^ ^ operators. The discretization results in one diagonal sub-block 
for each grid point for each of the three Each intermediate step of the algorithm 

is now premultiplied by the associated with the implicit operator. Writing the 

product operator ^ as L^, the inverse operators are thus removed and the scheme 

is written, once again dropping the vector superscript for convenience 


[q,a" - At/3L,"][4>»- $"] = AIL, 0"+ AtQ|[ij''+^"]4>" + AIQ,s"+^ 
[QjA"-A1/9Lj"][4>**-<I."] = QjA"[<I>‘- *"] 

[q,a"- AtaL,"][4.”*- <!."] = Q3 a"[4.**- *"] 




(j)n+i_ 0( At^) 
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With the removal of the inverse operator Q"^, the question of the proper inter- 
mediate level solution boundary conditions can be addressed. As is pointed out by 
Briley and McDonald (Ref. 2), the proper intermediate level boundary conditions may be 
derived by running through the intermediate steps in reverse order. Defining a 
boundary condition operator b” after linearizing the appropriate physical boundary 
condition by Taylor series expansion in time as 

and applying this operator to the algorithm defines the boundary conditions as 

=[b3"q3a'^ - At /3B3"L3"] 

(59) 

= [b2"Q2A" - 


Note that unless commute (an unlikely event except with Dirchilet boundary 

conditions, where B^ = I) the exact boundary conditions cannot be derived. A number 
of possible strategies are possible at this point aimed at various levels of approxi- 
mation to B^L^. For the present, the term AtgB^^[4> - is neglected. This 
introduces an error of order 0[At($ - $’^) ] into the solution but note that this error 
disappears at steady state where Neglect of the AtBB^^ [0 - term 

is, of course, equivalent to applying the physical boundary conditions on the inter- 
mediate level variables. 

This completes the general derivation of the algorithm and attention is now 
given to the specific forms of the operators including the rather special form 
of the component operator for the continuity equation. 

It is worth noting that the operator ^ or ..^can be split into any number of 
components which need not be associated with a particular coordinate direction. As 
pointed out by Douglas and Gunn (Ref. 19), the criterion for identifying sub-operators 
is that the associated matrices be "easily solved" (i.e., narrow-banded). Thus, 
mixed derivatives and the complicating terms which might inhibit the use of OCI can be 
treated implicitly within such a framework, although this would increase the number 
of intermediate steps and thereby complicate the solution procedure. 

An inspection of Eq. (58) reveals that only the linearized operators l”, and 
appear. Indeed, the computer code employs this feature by evaluating these three 


operators before the first sweep, storing them and accessing them as needed in the 
subsequent three sweeps. In addition, the terms arising from the nonlinear terms are 
immediately absorbed into as they appear, allowing for an efficient evaluation of 

the terms in the differential equations. 

The spatial operators appearing in the differential equations are 
2 must be identified at least formally in order to isolate the coefficients that 
are to be used in the construction of the Q and R operators. The operators 
^ 1* 2’ *^3 represented 'in standard form at each grid point, i.e., 


A " = a" <I) 


+a"4)+a"4> + a”tl> 

12 1,1 13 I 14 2 15 3 


(60) 


In Eq. (60) the first subscript of $ Indicates the velocity component (associated 
with the corresponding direction and " , " indicates a derivative. The subscripts of 
the a^j refer to the direction (1) and the term in the equation (j) respectively. 

Note that the equation is in quasi-linear form, since the coefficients of the 
derivative operators need to be identified, for use with the QR operator technique 
employed here. Alternate schemes have been proposed by Leventhal (Ref. 28) for 
equations in conservation form but are not considered here. In the following section, 
a description will be given of how this entire operator is discretized by employing 
the QR operator format, and how the discretization is incorporated into the LBI 
framework in order to solve the system of equations (58). 

The continuity equation is considered first. Since it is a first-order partial 
differential equation it does not have the standard form of Eq. (59). Furthermore, 
in the linearization process p has been eliminated in favor of the u^ velocity 
components so that the continuity equation has become an equation for the three 
velocity components, and not density. 

An inspection of the system of equations under consideration reveals that 
substantial savings can be realized if the equations are partioned appropriately. 

This is in keeping with the observations of McDonald and Briley (Ref. 8) who noted 
that skillful partioning of the resulting matrix can lead to significant decreases in 
computation time. Due to the use of a boundary layer coordinate system, the normal 
velocity appears only in conjunction with terms associated with the normal "3" 
direction in the two momentum equations. Hence, in the first two sweeps one is 
required to solve only for the two corresponding velocity components the streamwise 
and spanwise momentum equation without the need of considering the continuity equation. 
However, on the third sweep where all 3 velocity components appear, one must solve 
all 3 equations. This strategy reduces the solution procedure to the inversion of 
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two 2 X 2 block matrices and one 3 x 3 block matrix rather than three 3x3 block 
matrices which leads to a substantial reduction in computation time. If the full 
Navier-Stokes equations were considered (including a normal momentum equation) 
the aforementioned partloning could not be applied since the normal velocity would 
appear in all three sweeps. 

The question that arises is how to appropriately split the continuity equation, 
since it is solved only on the third sweep. Here again the Douglas-Gunn formulation 
leads to the appropriate choice. The continuity equation written in conservation 
form is, 



After linearizing and eliminating p, the increment form is obtained 


a"Au 


n+i 


At 

J 


+ b"Aw"'+ ^[vVau"*' + v"b'’Aw"*' + /Av"*'] 

■^[ ^ -^ [(/ +u"a'’)Au"*'+(uV)Aw"*'] 

+ ^ [(^" + „"b")Aw"*'+(wV)Au"*'] 


(62) 


where all the veloicty components are the contravariant components u = u^, w = u 
and V = u^. J. is the Jacobian end 


[a,,"" 9,2""] 

B" = •^[9aew''+g,j,u"] 


By employing the Douglas~Gunn procedure. Equation (62) is represented as a third 
sweep equation, and a consistent approximation is obtained to the continuity 
equation, i.e. , the x derivative term is evaluated at the * level and the x^ derivative 
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terra is evaluated at the ** level, The values of the Interroediate derivative terms 
are obtained after the solution of the first two sweeps of the two momentum equations. 
Note that these terms do not contain the normal velocity, The equation can thus be 
written in s;)nnbolic form 

J dx^ L i 


3 

Since the only term Involving v is in the x derivative term, one can directly 

3 

integrate the equation with respect to x , i.e. 



The next section describes how this is done very easily via the QR operator scheme. 

The concept of integrating directly the continuity equation is not now. Davis (Ref. 29) 
in his coupled procedure for the solution of two-dimensional steady boundary layer 
equation used a trapezoidal rule to integrate the continuity equation. Weinberg 
(Refs. 30 and 31) also used a fourth order Simpson integration scheme to solve the 
compressible boundary layer equations. Such procedures are stable and offer a viable 
alternative to approximatxng the derivatives by finite differences. Note that 
conceptually the continuity equation in integrated form is treated on each sweep 
of the Douglas-Gunn splitting, although in actuality this can be viewed as having 
the same form as each sweep and the integration operator can be incorporated into 
the and ^ difference operators, and as a result the stability and consistency 
of the original splitting is retained. 

Implementation of the LBI Scheme Employing 
the QR Operator Technique 

Consider the third sweep of Equation (57) in which both momentum equations and 
the continuity equation are solved. The momentum equations are in the form 
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ORfGiMAIL PAGE (9 
OF POOR QUALITY 


[a" A"Act** 

(64) 

where is the column vector of unknowns, u, v, w. Here it has been implicitly 

assumed that the equations have been appropriately normalized and that the contra- 
variant velocity components have been suitably transfonned into their physical 
components. Employing physical components, (cf. Ref. 32) leads to a better behaved 
solution since these components are not unduly influenced by geometrical variations. 
For the streamwise momentum equation one obtains 

Au, 33+ Ug3 Au,3 Q33Au + 0^3 Aw + OgjAv (65a) 

while for the spanwise momentum equation one obtains 

b23Aw,3 + b33Aw + b^sAv + bgjAu (65b) 

where superscript *** has been omitted from Au, Av and Aw. Now in Equation (65a), 
the first three terras on the right-hand side are approximated by the operator 
equivalent 


so that 



°r' 

Axj^ 


Au 




Q,”'r,Au 

— + a^aAw 4- 053 A V 


( 66 ) 


(67) 


Similar approximations are made for Equation (65b) . After substituting Equation (66) 
into Equation (64), and multiplying thru by Q, one obtains for the streamwise 
momentum equation 

[q,p" t )9XR|]Au -■ /3AI Q|0413Aw -/SA lOjOgjAv = QiP^Au**^ (68) 

2 

where X = At/Ax^ 

Similarly for the spanwise momentum equation, one obtains 

[Og/)^ - /3XRg] Aw - /SAtOgb^^gAv - ^SAtC^sjAu = Og^^Aw** (69) 

The same type of procedure is employed for the continuity equation. Since the 
continuity equation involves only first derivatives, they can be represented as 

a _ q~c'rc 

ax* " Axj 
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ORIGINAL PAGE 13 
OF POOR QUALITY 


The operotore Qc nml Rc ore constructed to npproximotc the weighto assoclnted with 
either n second order tropezoidal rule or o fourth order Simpson's rule, i.e, 


Trapezoidal rule 


Simpson's rule 




2 




1 




The discretized continuity equation thus becomes 

JA'’Au + Jb"Aw+ ■|^Q~'Rj.[dA%"Au+ JB%"A w + J/5"Av] = RHS (71) 

where RHS contains all the terras due to the linearl 2 ation procedure and the terras 

evaluated at the * and ** levels. Multiplying thru by Q and setting w * At/Axo, 

c o 

the equation reduces to 

[OgJA"'+ /3ujRc’dA"v"]Au + [Qj.JB"+ /3(uRgJB'’v^]Aw d*[ /ScuR^/j^JaV = Qg(RHS) (72) 

The resulting matrix derived from Equations (67), (68) and (71) becomes a block 
3 tridiagonal matrix (Q and R are tridiagonal operators) with each sub block taking 
on the form 


[ Q,/ -|GXR, ][ -/3AtO,a43 ][-^AlQ,a53] 


Au 


Qi(Au**) 

[ -/3AlQ2b53 ][ ][-^AtQ2b,J 


Aw 

= 

Q2(Aw*^ 

[OpJA^ +y9to R^JAV] [Qj,JB"+y6fcuRgJB%’] ] 


Av 


Q^(RHS) 

w 


This matrix is inverted by standard LU decomposition. 

Boundary Conditions and Initial Conditions 

The type of boundary conditions employed in the solution of the approximate form 
of the Navier-Stokes equations are described in this section. On the body surface no- 
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slip Is prescribed for all the velocity components. At the outer edge of the viscous 
layer the magnitudes of the streamwlae and spanwise velocity components are also 
prescribed, However, the value of the normal velocity component is not set, but 
rather computed as part of the numerical solution as is the practice in standard 
boundary layer procedures. 

At the inflow boundaries, (upstream) velocity prc files are fixed, while extra- 
polation conditions are employed at the outflow boundaries (dowstream) , Further 
discussion of this matter is given in the section of numerical results. 

The intermediate boundary conditions employed on the first two sweeps are the 
physical ones. For steady multidimensional problems, the imposition of physical 
intermediate boundary conditions did not impair the quality of the solutions 
obtained. TIaese results are in keeping with the analysis of McDonald and Briley 
(Ref, 2) for second order spatial schemes. For the unsteady cases considered 
physical intermediate boundary conditions have also been used without any apparent 
difficulty. 

The question of proper intermediate boundary conditions for fourth order 
methods until recently has not been resolved, Fairweather and Mitchell (Ref, 33 ) 
developed nonphysical intermediate boundary conditions for a fourth order solution 
of Laplace's equation, and showed that, in general, the use of noncorrected, i.e,, 
physical boundary conditions leads to a loss in steady state accuracy for their 
method. As pointed out by Fairweather and Mitchell (Ref. 33), their scheme is 
inconsistent. It is this Inconsistency that requires one to use appropriately derived 
intermediate boundary conditions in order to recover a steady state solution indepen- 
dent of time. However, if a consistent scheme were to be used, e.g. Douglas-Gunn, 
then physical boundary conditions can be applied without any loss in steady state 
accuracy, including fourth order generalized OCI schemes. These conclusions 
generalize the results obtained by Briley and McDonald (Ref, 2) for second order 
finite difference methods to higher order schemes and to those schemes that can be 
cast into a QR operator framework. However, using physical intermediate boundary 
conditions is expected to decrease the overall temporal accuracy to first order. 

The Computer Code 

The type of numerical algorithm employed as well as its formulation has a 
marked impact on the structure of the computer code. One needs to consider both 
the number of CPU operations as well as the memory requirements. Usually, the 
number of operations can be reduced at the expense of increasing the amount of 
storage. However, for three-dimensional problems the accessible fast (small core) 


memory becomes a severe limitation oven in the case where the code has not bet*n 
modified to optimize the operation count. 

the storage requirements for the solution of the approximate form of the three- 
dimensional Navler-Stokes equations for even modest size grids (e.g, 30 x 30 x 30) 
exceed the available small core memory of a machine like the GDC 7600 r One must then 
resort to either mass storage devices such as disks or slow access tceraory (large core) . 

In using such devices both access time and transfer rates must be considered. 

When small amounts of data are being transferred frequently then access time becomes a 
significant factor. Therefore, a combination of strategies must be employed in order 
to optimize both access time and transfer rate. 

An investigation of the operation count of the LBI shceme in conjunction with 
the QU operator technique leads to the conclusion that the most significant fraction 
of time is spent in computing the maLrix coefficients, l.e. the linearization coeffi- 
cients and difference weights. This amount far exceeds the time required for the matrix 
inversion. Hence, it is worthwhile to optimize the calculation of these coefficients, 
and if possible store their values. This procedure was accomplished by storing the 
operator coefficients 2 and ^ ^ they were computed in the first sweep on the 
right-hand side of the differential equation. On the second and third sweeps ,^ 2 
and o were accessed respectively and were not recomputed. It was for this reason 
that the formulation of the LBI scheme referred to the linearized operators ^ 

Instead of the t> ^'s on the right-hand side of the equation. 

In order to minimize data transfer, the code was constructed to have in memory 
one plane of data at any given time. Recall that 1-2 planes are parallel to the 
surface while the 3 direction is normal to the surface. During the solution procedure, 
which is described below, the first and second sweeps are conducted on 1-2 planes, 
evaluated in sequence in the 3 direction from the surface to the outer edge. The 
third sweep is conducted on 1-3 planes evaluated in sequence in the 2 direction, 

A problem arises in two dimensions during the first sweep, in that only the '1' 
or streamwise direction is required yet '1-2' planes must be solved sequ.i iitially. This 
means that only '1' lines are solved in the first sv/eep since the '2' direction is 
essentially passive, and thus leads to unnecessary data transfer. In order to 
alleviate this inefficiency the code was modified so that in two dimensions the 
equations would be solved 'in core' in the '1-3' plane on both sweeps. At the same 
time, additional modifications were introduced that resulted in a significant speed 
up of the code. These include specialized hard wired 2x2 and 3x3 block matrix 
multipliers and inverters, and an efficient addressing routine for spatially 
deper’o.nt variables. As a result of Including these modifications, the CPU time for 


- m 




two-dimansionaX problems was reduced by nearly an order of magnitude, The run time 
for a two-dimensional turbulent calculation is ,00063 sec per grid point per time step 
and for a laminar two-dimensional calculation is ,00035 sec per grid point per time step. 
Although there are still additional modifications remaining to be incorporated into 
the code, it is believed that the major gains in efficiency have already been obtained. 

The general structure of the computer code will now be described. After the 
input section and the initialization of data e,g. geometry, grid transformations, 
flowfield,. etc. the actual construction of the difference operators is begun. The 
first derivatives of the velocity components and viscosity are obtained for the entire 
flow field and stored for ready access when needed for the computation of the appro- 
priate terms in the governing ec^’ations. Tliereafter the terms that are to be treated 
explicitly are evaluateu and abu rbed into the function . 

The operators y "^2 *^3 then computed. These are used to 

evaluate the appropriate Q and R coefficients which are then stored for easy retrieval 
during each of the ADI sweeps. 

In the first sweep the matrix resulting from the application of the 
operators for the streamwise and spanwise momentum equation is solved as a 2 x 2 
coupled system. The solution of this system, the * level quantities, are then used 
to construct the right-hand side of the second sweep equations and to evaluate the 
appropriate * level term in the continuity equation. At this point the ^2 op^^^^tor 
is accessed and again a 2 x 2 system of equations for the streamwise and spanwise 
momentum equation is solved. The ** level quantities are then used to construct the 
right-hand side of the third sweep equations as well as the appropriate terras in the 
continuity equation. For the third sweep equations which consist of the two momentum 
equations and the continuity equation, the 0^2 oP®^^tor is accessed from memory. The 
resulting 3x3 system of equations is solved for the three velocity components. 

After the primary variables are evaluated, the theromodynamic quantities, density, 
temperature and viscosity are computed. The procedure is then repeated at the following 
time steps. 

It is noteworthy that the scheme just described operates on vectors, i.e. lines 
of data. Therefore, it could show promise for vectorized machines. 
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Numerical Results 


In this section we describe results of numerical computations that were ob- 
tained by exercising the computer cotte discussed previously. The objectives of 
these calculations are; i demonstrate the viability of the code in performing two 
dimensional steady turbulent flow calculations as well as in computing two dimensional 
unsteady flows, to compare second order finite differences with fourth order methods, 
to discover any limitations of the solution procedure and to indicate what modifications 
would be recommended to improve the code's performance. 

In meeting these goals initial calculations were performed with the fourth order 
generalized OCI scheme, Model one and two dimensional problems were considered, 
progressing to the solution of the laminar boundary layer equations. These calcula- 
tions which included both the Blaslus and Howarth flows were used to validate the 
fourth order generalized OCI option as well as the code modifications described in 
the previous section. 

At this point the choice of the specific calculations to be considered was 
addressed. The first set of calculations consisted of two turbulent cases; zero 
pressure gradient (Weighardt Flat Plate (Ref. 34)) and adverse pressure gradient 
leading to separation (Newman Airfoil (Ref. 34)). These cases were chosen due to 
the reliable experimental data that are available with which comparisons can be made. 

The second set of calculations treat the laminar two dimensional oscillating flow 
over a flat plate. Here again there is abundant information (theoretical as well as 
computational) with which to make meaningful comparisons. 

Two-Dimensional Steady Turbulent Flow 

Weighardt Flat Plate 

The first case we consider is the turbulent flow over a plat plate with zero 
pressure gradient. Flow conditions were specified to match with the experimental data 
of Weighardt (Case 1400 in Ref. 34), viz. 

n = 108 ft/sec. 
e 

R^ = 661,190/ft. 

The computational domain was chosen with the inflov; boundary located at 
X = .6135 ft and the outflow boundary located at x = 13.1232 ft, while the outer 
edge was set to a constant value of y = .25 ft for all x. In order to resolve the 
sublayer two types of grid transformation were employed, a Roberts type (Ref. 11) 
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which is essentially a hyperbolic tangent function and an Oh type (Ref. 35) which 
is constructed from a series of complimentary error functions. In the streamwise 
direction a uniform grid distribution was employed. The total grid distribution 
consisted of 36 points in each direction. 

At the inflow boundary the velocity profiles are required to be specified. 
These were obtained directly from Weighardt's data by using a Coles' law profile 
(Ref. 34) l.e. 



+ C + 


27T(x) , 2 

— ^sin 



I Wall-Wake Law 



laminar sublayer 


where 


y* 


yur 

u 



Ur " 


(73) 


and 7t(x)/k is evaluated from the condition that u = u^ at y/6 = 1. In Eq. (73) the 

constants k and C are set at .41 and 5.0 respectively. This leaves two free parameters, 

T (or C^) and 6, which must be chosen to completely specify the profile. In (Ref. 34) 
w t 

tabulated values of and 6 are given, which are evaluated from curve fits of the 
experimental data as a function of streamwise position. The values corresponding to 
X = .6135 ft are 

= .004138 
6 = .015472 ft. 

Once the streamwise velocity profiles are obtained the normal velocity can be 
determined by approximating 9u/8x and then integrating the continuity equation. 

In the calculation procedure the flow was assumed incompressible and in the 
streamwise direction the boundary layer option was employed, i.e. streamwise dif- 
fusion terras were dropped and a backward difference approximation was used for the 
streamwise convective terms. In the direction normal to the wall both second order 
finite differences and the fourth order G/OCI schemes were exercised. The full form of 
the dissipation function is used in the evaluation of the turbulent viscosity since 
in this case it did not have any perverse behavior. In the Newman airfoil calculation 
to be discussed shortly, however, the uSe of the full dissipation function adversely 
influenced the behavior of the solution. The solutions were obtained by marching in 
time until a steady state was reached. 


41 




# 



The second order and fourth order calculation required 42 and 43 Iterations 
respectively to converge, where a convergence criterion of e«10~^ was employed and e 
corresponds to the maximum change in a primary variable (u,v) over a range of time steps. 
The total run time for the fourth order G/OCI calculation was 35.28 sec on a CDC 7600 
computer which reduces to .00063sec per grid point per time step. 

The skin friction distribution, as a function of distance along the plate is 
presented in Figure 1. A comparison between the experimental data and both computa- 
tions shows good agreement over most of the plate, with the major discrepancy appearing 
near the upstream boundary. There are two reasons for this. First, there is an 
incompatibility between the inflow profiles (obtained indirectly from experimental 
data) and the velocity profiles downstream of the inflow plane, which are obtained 
by the numerical solution of the governing equations. There is no direct method of 
avoiding this inherent error other than refining the mesh in the streamwise and 
normal directions in the vicinity of the inflow plane. This would allow any errors 
that are generated at the boundary to damp out several grid points downstream. 

Another method that could be used relies on completing two computations; the first 
beginning upstream of the inflow plane and extending downstream of it, and the second, 
the one of interest, on the actual computational domain. In the second calculation 
for the inflow profiles, one uses the profiles computed from the first calculation at 
the corresponding x location. Although such a procedure is not always applicable, it 
was used successfully for the Newman Airfoil Case (described below) . 

The second cause of the discrepancy, and probably the more serious, was the lack 
of resolution of the velocity profile at the upstream plane. Recall, that the top 
boundary was prescribed at some fixed constant height for all streamwise locations. 

The extent of the normal distance is based upon the consideration that at the down- 
stream boundary the velocity profile should be totally contained within the computa- 
tional domain. As a result at the upstream boundary the "boundary layer" occupies a 
small portion of the entire normal extent, leading to a lack of resolution there. 

This problem can easily be remedied by introducing ay/6 transformation which accounts 
for boundary layer growth and assures equal resolution in the normal direction for 
all X locations. In view of the benefits of this type of transformation it is 
recommended that the computer code be modified to include a y/6 transformation 
on option. 
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Despite the shortcomings discussed above, the errors near the inflow boundary 
die out rapidly and for the remainder of the flow field computations are in excellent 
agreement with the data. The velocity profile at x « 11.44 ft. (of Figure 2), also 
shows that good agreement was obtained between the computations and experiment. 

In Figure 1, the increased accuracy of the fourth order generalized OCI scheme 
is evident. However, there were qualitative differences between the second order and 
fourth order solutions. The streamwlse velocity profiles obtained by second order 
finite differences were smooth throughout the entire boundary layer. In contrast to 
this behavior the fourth order generalized OCI solutions tended to have slight wiggles 
near the outer edge of the boundary layer even though in theory the generalized OCI 
scheme should not if certain inequality constraints are satisfied. Unfortunately, 
near the outer edge these conditions are violated allowing for the observed behavior. 

We believe that the wiggles are due in part to high order numerical differentiation 
of the turbulent viscosity which possesses a sharp knee in the wake region of the 
layer. It is well known that higher order derivatives of rapidly Varying functions 
introduce noise, and it is this noise which we believe we are witnessing. Laminar 
flow calculations have demonstrated that oscillatory behavior is not exhibited near the 
outer edge even for very coarse meshes. It is, therefore, felt that in order to obtain 
turbulent solutions commensurate with laminar solutions additional investigations of 
'smoothing' procedures for the turbulent viscosity should be undertaken, and appropriate 
grid distributions be considered. 

As a final note, the computation of C^ is discussed. Second order finite dif- 
ferences employ a three-point second-order one sided difference. A comparable formula 
accurate to fourth-order requires a five-point one-sided difference. As a result of 
the grid transformation employed, the formula encompasses points relatively far from 
the wall, and hence the computed will be in error. Instead, for the reported 
results In Figure 1, a three-point, one-sided formula was used which gives more accurate 
values than the five-point scheme. Such behavior is not uncommon in computing one- 
sided derivatives (cf. Ref. 36). 
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Newman Airfoil Case 

The Newman airfoil flow (case 3500 - series 2, Ref, 34) is considered next. It is 
more interesting and more difficult than the flat plate case due to the adverse pres- 
sure gradient present which leads to flow separation near the trailing edge. Indeed, 
the only meaningful numerical results obtained were those computed by second-order 
finite differences. The fourth-order calculations computed on the same grid were 
inaccurate and exhibited excessive oscillations near the outer edge of the computational 
domain. An explanation is suggested for the observed behavior, and a discussion of the 
computed results is presented. 

In the computations a 36 x 36 grid was employed. In the streamwise direction 
which extends from x ® 2.958 ft. to x = 4.926 ft. the spacing was uniform, while in the 
normal direction which extends from the wall to a fixed outer edge at y = .35 ft., 
a Roberts type stretching was used. The stretching parameter was chosen so that at 
the inflow plane the first grid point would be located within the sublayer (y = 1.75). 
The Reynolds* number for the flow was 769,231 per foot. 

Initial flow computations revealed that the results were sensitive to both the 
external velocity distribution and the upstream velocity profile. Differentiating the 
velocity data given in Ref. 34 directly for use in the computation of the pressure 
gradient was found to work best. Numerically differentiating analytical curve fits 
to the external velocity distribution Introduced errors in the pressure gradient and 
were thus discarded. In order to eliminate incorapatabilitles caused by a mismatch in 
the inflow conditions and the numerical solution, a preliminary calculation was per- 
formed on a 21 X 36 grid in the domain 2.759 < x <3.2, and the computed profiles at 
X = 2.958 were used as upstream conditions for the primary calculation. Since the 
pressure gradient in the upstream region was mild, there was no difficulty in obtaining 
a converged calculation for the 21 x 36 grid (the calculation converged in 21 iterations). 
However, the "boundary layer" form of the dissipation function was used since with that 
form the solutions behaved better near the outer edge of the computational domain. 

In Figure 3 the computed skin friction distribution is presented and compared to 
the experlmentivi data. The agreement is good for most of the airfoil with the exception 
of the trailing edge region where the second-order calculation predicts separation up- 
stream of the actual separation point. However, for the shape factor, H = 6 /0 (cf. Fig. 4) 
the agreement between data and computations is not as good, with the numerical results 

X 

underpredicting the data. A lower value of H signifies that 5 is too small, 6 is too 
large or a combination of the two. In any event, the lack of agreement in H means that 
the velocity profile shape is incorrect. Two effects were considered to account for 
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this discrepancy; lack of spatial resolution and turbulence modeling. As discussed in 
the previous section, the inability to account for boundary layer growth leads to a 
spatial resolution problem. This effect manifests itself particularly near the upstream 
boundary where the boundary layer is thinnest. Since the major discrepencies in this 
case lie near the downstream boundary, where the boundary layer is thickest, poor grid 
resolution caused by the lack of a y/6 transformation is not a predominant factor for 
the observed behavior. Of course, grid resolution could be a real concern, but only in 
how it affects the numerical scheme in resolving regions of large gradients, l.e. the 
knee in the viscosity profile. In order to investigate the effect of turbulence model- 
ing, a series of calculations were performed using a modified version of Cebeci's 
turbulent boundary layer code (Ref. 37) which uses Keller's Box scheme. This code per- 
mitted us not only to investigate turbulence modeling effects, but afforded us with a 
means of comparing our computational procedure with a different method, and running a 
series of mesh refinement studies. 

The modified Cebeci code was checked out by computing the Weighardt case. After 
successful completion of this calculation, the Newman Airfoil case was attempted. In 
running the Newman Airfoil anamolies appeared which were not present in the Weighardt 
case. The computed skin friction coefficient exhibited oscillations in the streamwise 
direction. Since the "mean" curve compared well with our calculations for x > 3.5 ft. 
and since our Intent was to investigate the effects of turbulence modeling and grid re- 
finement, we did not pursue the anamolous behavior. Even though oscillation in were 
observed, we believe that they do not invalidate the conclusions to be drawn. The 
calculations were performed with 36 points in the streamwise direction (matching our 
calculations), and either 36, 60 or 100 points in the normal direction, which was 
stretched by a logarithmic transformation. It is important to point out that these 
calculations were performed in similarity variables so that the growth of the boundary 
layer was taken into account. 

The results of the Box Scheme calculations were rather surprising. For the 

calculations using either Cebeci's two layer mixing length model (Ref. 14) 

or McDonald's model, the distribution varied insignificantly from our second-order 

finite difference solution shown in Figure 3. For the comparison, the oscillating 

part was neglected and a mean curve was considered. Neither was there much 

Improvement with the use of a finer mesh indicating that 36 point is probably sufficient 

to resolve the flow. The major differences between Cebeci's model and Equation 37 

* 

are that in Cebeci's model, the outer layer mixing length is proportional to 6 , while 
in Equation 37 it is proportional to 6 , and in the Van Driest damping Cebeci includes 
a t-erm that accounts for the pressure gradient. 
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Tha effects of these differences were considered next. In order to investigate 
the effect of pressure gradient terms in the Van Driest damping, the correction term 
was Included in McDonald's model (Eq. 37). The results of the calculation with 36 
points in the normal direction revealed that there was a significant shift of the 
curve upward coward the data, with much better agreement obtained with the data, 

However, the oscillation still persisted. Why that term did not have a similar effect 
in Cebeci's mixing length model we cannot answer, and therefore, we believe further 
investigation is warranted. 

In comparing the shape factor, the H distribution exhibited different behavior, 
These results are shown in Figure 4 which presents H as a function of x. It is apparent 
that the calculations using Cebeci's mixing length give better agreement with the data. 
Since the skin friction distribution did not compare as well, the good agreement of the 
shape factor with data must be related to the overall shape of the velocity profile. 
Hence, the velocity profiles obtained from Cebeci's model and from Eq. 37 were compared 
with data at x = 4.509 ft., (cf Fig. 5) to discover how the turbulence model effects 
the shape and thus H. Although the profile using Cebeci’s model fits the data more 
closely, the shape of the velocity profile appears to be in error. Therefore, the 
good agreement for H shown in Figure 4 may be fortuitous. Further investigations 
beyond the scope of the present effort would be necessary to draw additional conclusions. 

The following conclusions can be drawn from the calculations employing the Box 
Scheme: 

(a) Mesh resolution in the normal direction was not a significant factor in 
the observed results, 

(b) Choice of turbulence model can have a significant effect on the solutions. 

In view of the above, the inclusion of a one-equation and/or a two-equation turbulence 
model in the computer code is recommended in any future effort. 

Until now, we have neglected to say anything about the fourth-order calculations. 
The results obtained for a 36 x 36 grid exaggerated the mild oscillations in the 
velocity profiles observed in the Weighardt case. The inequality constraints required 
for nonoscillatory behavior of the generalized OCI scheme were violated. In order to 
satisfy these constraint conditions, more grid points, better transformations or both 
are necessary. It is, therefore, recommended that further investigation be performed 
in this area. 

Two-Dimensional Unsteady Flow 

The unsteady flow case considered is that of a nonzero mean flow with a sinusoidal 
unsteady component superimposed upon it, i.e. 
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Ug(t) = UqII +ACOSWt) 


(74) 


where is the mean flow velocity, A is the dimensionless amplitude of the periodic 
part and lo is the circular frequency. For the case of laminar flow over a flat plate, 
Llfththill (Ref. 38) obtained expressions for the skin friction coefficient in the limit 
of low and high reduced frequency. 

(0,332 + A(0,498COSojt - 0.849o)X/Uo sin lut ) OJX/Uq «l 

(75) 

0.332 -i a(u)X/Uq)'^^ COS(tot + tt/4) ojX/Uq » I 

It is convenient in the flow analysis to calculate the phase angle between the 
external velocity and a boundary layer property (l.e, skin friction coefficient) as a 
function of streamwise location, x « x^. The derivation of the phase angle is 
presented below. 

Consider an unsteady flow with an external velocity field that is a function of 
X and t 

Ug(x , t) = Uq(x) ( I +A COSojt) (75) 

and which has associated with it a boundary layer property f(x,t). Denoting the 
average value of a function by a super bar (~) , the average external velocity at 
point x^ is defined as 

I I r i r i 

u(Xo)= — J Uo(xo,t)dt [ Uo(Xo)J [1 + a cos aj?J dt (77) 

or 

il(Xo) = Uq{Xo) (78) 

where T = 27r/o) is the period of oscillation ant t^ is a reference time from which 
averaging begins. Similarly, the average value of the function f becomes 

f(Xo) = Y /*'^Hxo,t)dt (79) 
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Here we have tactly aasumed that the flow has reached a pseudo-steady state and 
initial transients have died out, For the cases considered here, such a state is 
reached after one period so that tj^ » T, 

The function f(x,t) is now expressed as a periodic function with frequency u, 

i.e, 

f (x,t) = 7(Xq) + Bcos[(ot + ( 80 ) 


where is the phase angle between u^ and f at point x ■ x^. 

In order to determine the phase angle, the following Integrals are evaluated 


cos 


[(^(Xq)] 


A. 

T 


J|+T 


J\i ’‘•^o(Xq)^« (f{XQ,t) - fCXq)) dt 

B AUo(Xq) 


and 


B 



/'(f(x„,l) -T(Xj,)) 


2 dt 


(81) 


(82) 


The character of the flow field is governed in part by the amplitude of oscil- 
lation, A. Ackerberg and Phillips (Ref, 39) have shown that for A sufficiently small 
A < .3 no backflow will occur. However, for larger values, A > ,3 regions of reverse 
flow will appear over the plate, predominantly near the downstream boundary. The 
occurrence of reverse flow at the downstream boundary has serious implications on the 
numerical solution of the governing equations, and will be discussed in greater 
detail below. 

For the case considered in this report, A was chosen so as to avoid reverse 
flow from occurring, while the other input quantities were specified by numerical 
considerations (see below). Hence, the input data employed in the calculation are 


u = .10 
o 

A = .125 
u) = 5n/2 

Xy « .1 

^ = 3.6 
Re = 666,667 
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. .V 


Before describing the parclculars of the colculatlons, a general 'Uscusaion 
of the solution procedure, Its features, and its associated constraints will be given. 
The solution procedure for the unsteady case follows along the same lines as for the 
steady case, with the essential difference being that in the unsteady case, the 
solution is advanced in real time with the time step chosen from considerations of 
temporal accuracy rather than rate of convergence. Whereas in the steady case, all 
boundary conditions are time invariant, for the oscillating flow the velocity distribu- 
tions at the upstream and outer edge boundaries change with time. TIk is no dif- 
ficulty in applying the outer edge velocity boundary condition since the velocity 
distribution is simply that which is given in Eq. (74). The upstream boundary is 
more troublesome. If the upstream boundary is located at x * 0, then at that point 
the boundary layer is of zero thickness and hence in transformed boundary layer 
coordinates remains fixed. In that case, the upstream inflow boundary condition remains 
Blasius independent of time. However, if the upstream boundary is located at some 
small, but finite, distance downstream of the leading edge, the velocity profile will 
change with time. Since the calculations were conducted in the physical plane rather 
than in the transformed plane, the upstream boundary was required to be placed at a 
small finite value of x. Hence, in order to account for the varying upstream boundary 
layer profile a method used by Singleton and Nash (Ref. 40) was employed, viz. the 
upstream boundary was scaled by the new edge velocity at each new time level. This 
procedure fixes the values of the flow variables at the upstream boundary, and permits 
the use of function conditions (necessary for well posedneas) but as a result also 
introduces errors there. 

By solving in the physical plane, boundary layer growth could not be very satis- 
factorily accounted for. This resulted in having the outer edge fixed at y = 0.037 ft.; 
the choice of this distance being predicated upon the necessity to accommodate the 
boundary layer at the downstream boundary. In Figure 6 the computed phase angle <}) is 
presented and compared with Cebeci's calculations (Ref. 41), and to the low and high 
frequency predictions of Llghthill (Ref. 38). Both second-order and fourth-order 
calculations are shown, for 0=1. The agreement of our computations with the other 
predictions is very good for m _> .8. For small x there is a discrepancy due to the 
implementation of the upstream boundary condition described above, and the lack of a 
y/6 transformation. This effect was studied in greater detail at a lower frequency, 

0 ) = tt/ 2 and were compared to the calculations of Murphy (Ref. 42). 
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In Figure 7 second-order solutions with g ■ 1 and (.i » 1/2 (Cronk Nieolson) are 
compored to Cebcei's results. The Crank Nieolson calculation gives larger values 
of (f> than both the results for 6*1 and Ref. (41). Since Murphy's and Prenter's 
higher order calculations (Ref. 43) show the same trend, this would appear to 
indicate the improved accuracy of the Crank Nieolson scheme. 

As A is Increased, Phillips and Ackerberg (Ref. 39) show that reverse flow 
will occur during the cycle over some portion of the flat plate. Including the 
downstream boundary. Although unsteady boundary layers permit regions of reverse 
flow, there is a difficulty in applying appropriate boundary conditions if at the 
downstream boundary fluid is entering the computational domain rather than exiting it. 
In an attempt to remedy this situation, Ackerberg and Phillips argue that since 
disturbances travel at some finite speed then at a distance sufficiently far down- 
stream the flow will not have felt the disturbances generated at the leading edge 
and would, therefore, appear to be that which would exist on an infinite flat plate. 
The solution to that problem is well known, i.e. Rayleigh flow. Hence, they suggest 
that the Rayleigh solution be set at the downstream boundary (function condition) 
consistent with the corresponding x and t. 

We investigated the behavior of our numerical scheme for a problem where reverse 
flovj- occurs, A = .3. The boundary conditions described in Ref, (39) were used for 
the solution of the full equations and central differences with artificial dissipa- 
tion was fcmployed in the streamwise direction. The calculations were run for two 
cycles, but the transients had not as yet died out completely, The results which 
we are not presenting at this time (since they are of a preliminary nature) Indicate 
that there was no difficulty in the computations even when there were large regions 
of reverse flow. Future work will be aimed at conducting a more comprehensive study 
in this very important area. Additional effort will be extended on turbulent 
unsteady flows. 
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CONCLUSIONS 


In this report we have presented results obtained by exploiting the computer code 
developed under a previous phase of the research effort to solve an approximate form of 
the time-dependent Navler-Stokes equations, Both two-dimensional turbulent and unsteady 
laminar cases are considered. The governing equations that are solved are more general 
than the conventional boundary layer equations, notably in the inclusion of streamwlse 
and spanwlse diffusion terms, although the pressure is still imposed by the external 
flow, os in conventional boundary layer theory. The computer code incorporates the split 
LBI scheme in conjunction with OR operator scheme that permits a variety of spatial dif- 
ference schemes, including standard second-order finite differences, exponential type 
methods and fourth-order OCI techniques. In the split LBI scheme, an implicit sweep is 
performed in each spatia] coordinate direction, A careful ordering of these sweeps 
permits an uncoupling of the continuity equation from the system in the first two 
implicit sweeps. Thus, on the first two sweeps the (trldlagonal) system block size is 
reduced from 3 x 3 to 2 x 2 with a resulting cost savings. On the last sweep of each 
time step all the equations in the system are linearly coupled and 3x3 blocks must be 
eliminated . 

Results of computation Indicate that the procedure is viable far more complex 
problems of Interest. Higher order methods can yield more accurate results although 
care must be taken fo? turbulent flows when coarse grids are employed. For unsteady 
flows, the method is extremely efficient as a result of the noniterative nature of the 
algorithm. Future efforts will be aimed at incorporating ay/6 transformation, which 
will eliminate some of the difficulties encountered and will make the code more robust. 
Further investigations of appropriate turbulence models is also recommended with the 
one equation x-a model given priority. 
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AI’PEf.’DlX A 


Linearization Technique 

A number of techniques have been used for implicit solution of the 
following first-ordev nonlinear scalar equation in one dependent variable 
•J-Cx.t): 


d0/dt =F{<p) dG(<^)/dx 


(Al) 


Special cases of Eq. (Al) include the conservation form if F((f') ■= 1, and 
quasi-linear flow if G((J)) = Previous implicit methods for Eq , (Al) 
which employ nonlinear difference equations and also methods based on two- 
step predictor-corrector schemes are discussed by Ames (Ref. 44, p. 82) and 
von Rosenburg (Ref. 45), p. 56). O.ne such method is to difference nonlinear 
terms directly at the implicit time level to obtain nonlinear implicit 
difference equations; these are then solved iteratively by a procedure such 
as Newton's method. Although otherwise attractive, there may be difficulty 
with convergence in the iterative solution of the nonlinear difference 
equations, and some efficiency is sacrificed by the need for iteration. An 
implicit predictor-corrector technique has been devised by Douglas and Jones 
(Ref. 46) which is applicable to the quasilinear case (G = ij>) of Eq . (Al) . 

The first step of their procedure is to lineariJte the equation by evaluating 
the nonlinear coefficient as F(i}i”) and to predict values of using either 

the backward difference or the Crank-Nicolson scheme. Values for ()> are 
then computed in a similar manner using F((fi’^'^^^^) and the Crank-Nicolson scheme. 
Gourlay and Morris (Ref, 47 ) have also proposed implicit predictor-corrector 
techniques which can be applied to Eq. (Al) . In the conservative case (F = 1) , 
their technique is to define G(<{)) by the relation G((j)) = <^G((j)) when such a 
definition exists, and to evaluate G(i}) ) using values for <ji computed by 

an explicit predictor scheme. With G thereby known at the implicit time level, 
the equation can be treated as linear and corrected values of are computed 

by the Crank-Nicolson scheme. 

A technique is described here for deriving linear implicit difference 
approximations for nonlinear differential equations. The technique is based 


on an expansion of nonlinear implicic terms about the solution at the known 
time level, and leads to a one~step, two-level scheme which, being linear 
in unknown (implicit) quantities, can be solved efficiently without iteration. 
This idea was applied by Richtmyer and Morton (Ref. 15 > p* 203) to a scalar 
nonlinear diffusion equation. Here, the technique is developed for problems 
governed by 1 nonlinear equations in 1 dependent variables which are functions 
of time and space coordinates. The technique will be described for the three- 
dimensional, unsteady equations. 

The solution domain is discretized by grid poinl:.s having equal spacings 

12 3 12 3 

in the computational coordinates, Ay , Ay and Ay in the y , y and y 

directions, respectively, and an arbitrary time step. At. The subscripts i, j, 

12 3 

k and superscript n are grid point indices associated with y , y , y and t, 
respectively, and thus denotes <J>(y^, y^, y?, t”) . It is assumed that 

1,J,K 1 J K 

the solution is known at the n level, t^, and is desired at the (n+1) level, 
n*f*l 

t . At the risk of c.n occasional ambiguity, one or more of the subscripts 
is frequently omitted, so that is equivalent to . . . 

The numerical method employed is quite general and is formally derived for 
systems of governing equations which have the following form: 

■ (A2) 


where ({i is a column vector containing dependent variables, H and S are 
column vector functions of (p, and ^ is a column vector whose elements are 
spatial differential operators which may be multidimensional. The generality 
of Eq. (A2) allows the method to be developed concisely and permits various 
extensions and modifications (e.g., noncartesian coordinate systems, turbulence 
models) to be made more or less routinely. It should be emphasized, however, 
that the Jacobian must usually be nonsingular if the ADI techniques as 

applied to Eq. (A2) are to be valid. A necessary condition is that each 
dependent variable appear in one- or more of the governing equations as a time 
derivative. An exception would occur if for instance, a variable having no 
time derivative also appeared in only one equation, so that this equation could 
be decoupled from the remaining equations and solved £ posteriori by an alter- 
nate method. 
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The linoarlzod difference’ approximation is derived from tlit following 
implicit time-difference replacement of Eq. (A2) : 

where, for example, The form of and the spatial differ- 

encing are as yet unspecified. A parameter 6(0 £ 6 - 1) has been introduced 
so as to permit a variable centering of the scheme in time. Equation (A3) 
produces a backward difference formulation for 6 = 1 and a Crank'-Nicolson 
formulation for 6 = 1/2. 

The linearization is performed by a two-step process of expansion about 
the known time level t” and subsequent approximation of the quantity 
(3(}>/3t)^*At, which arises from chain rule differentiation, by - <})"). The 

result is 

= ( 9 b" + '-^") + 0(/it)* (A«a) 

. S'’+'=s"+(as/d^.)" (A«b) 

aw"’"') = ^(^")+(a^/a^W"'^'-'f") + o(Ai)" (A«c) 


The matrices 3H/3(fi and 3S/3(f’ are standard Jacobians whose elements are defined, 

for example, by (3H/3tf0^j. = 3H^/3(|>^. The operator elements of the matrix 

3^/3<|i are similarly ordered, i.e., (3^/3(j)) = 3.2) /9(j> ; however, the 

r cj ^ 

intended meaning of the operator elements requires some clarification- For 


the q'"^ row, the operation ( 3.2) 3(f)) ”((j-’^ - is understood to mean that 

{ 3/3t2> [(}i(x,y,z,u) ] }”At is computed and that all occurrences of (3(|)j-/3t)^ 


,n+l 


arising from chain rule differentiation are replaced by (cf.^ - 4>^)/At. 

After linearization as in Eqs. (A4) , Eq. (A3) becomes the following linear 
implicit time-differenced scheme: 
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n+1 

Although H is linearized to second order in Eq. (A^) , the division by At 
in Eq. (A3) introduces an error term of order At. A technique for maintaining 
formal second-order accuracy in the presence of nonlinear time derivatives is 
discussed by McDonald and Briley (Ref. 8), however, a three-level scheme 
results. Second-order temporal accuracy can also be obtained (for B “ 1/2) by 
a change in dependent variable to if ~ H(<J)) , provided this is convenient, since 
the nonlinear time derivative is then eliminated. The temporal accuracy 
is independent of the spatial accuracy. 

On examination, it can be seen that Eq. (A5) is linear in the quantity 

Ti 

(if) - (f* ) and that all other quantities are either known or evaluated at 
the n level. Computationally, it is convenient to solve Eq. (A5) for 
(ij) - (|) ) rather than . This both simplifies Eq. (Af) and reduces 
roundoff errors, since it is presumably better to compute a snail 0(At) change 
in an 0(1) quantity than the quantity itself. To simplify the notation, a 
new dependent variable \p defined by 

(A6) 

is introduced, and thus arid = 0. It is also convenient 

to revjrite Eq. (A5) in the following simplified form; 

(A+ Af = At [j^ (‘^")+s"] 


where the following symbols have been introduced to simplify the notation: 

A= -ySAl (3S^/c3(^) (A7b) 

J! ^ - fB{d :b /d<p) (A7c) 

It is noted that is a linear transformation and thus -^(0) = 0. Further- 
more if is linear, then J’(i|>) = -6^(i|/). 
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Spatial differencing of Eq. (A7a) is accoraplished simply by replacing 

i 2 i i 

derivative operators such as 3/3y , S /Jy 9y by corresponding finite 

2 fi 

difference operators, D^, D^. Henceforth, it is assumed that and J. have 
been discretized in this manner, unless otherwise noted. 

Before proceeding, some general observations seem appropriate. The 
foregoing linearization technique assumes only Taylor expandability, a'n assump- 
tion already implicit in the use of a finite difference method. The governing 
equations and boundary conditions are addressed directly as a system of coupled 
nonlinear equations which collectively determine the solution. The approach 
thus seems more natural than that of making £d hoc linearization and decoupling 
approximations, as is often done in applying implicit schemes to coupled 
and/or nonlinear partial differential equations. With the present approach, 
it is not necessary to associate each governing equation and boundary condition 
with a particular dependent variable and then to identify various "nonlinear 
coefficients" and "coupling terras" which must then be treated by lagging, 
predictor-corrector techniques, or iteration. The Taylor expansion procedure 
is analogous to that used in the generalized Newton-Raphson or quasi- 
linearization methods for iterative solution of nonlinear systems by expansion 
about a known current guess at the solution (e.g., Bellman & Kalaba, Ref. 48). 
However, the concept of expanding about the previous time level apparently 
had not been employed to produce a noniterative implicit time-dependent scheme 
for coupled equations, wherein nonlinear terms are approximated to a level of 
accuracy commensurate with that of the time differencing. The linearization 
technique also permits the implicit treatment of coupled nonlinear boundary 
conditions, such as stagnation pressure and enthalpy at subsonic inlet 
boundaries, and in practice, this latter feature was found to be crucial to 
the stability of the overall method (Ref. 49). 
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APPENDIX B 


GEOMETRIC PROPERTIES OF AIRFOILS 


original page 13 
OF POOR quality 


In this section we consider the vector properties of curves lying on a surface 
and relate them to their tonsor equivalents. In particular a NACA four-digit airfoil 
section is considered, and its geometric properties are computed in vector and tensor 
form, 

BASIC CONCEPTS 

Consider a radius vector r drawn from an inertial reference frame 0 to a point 
P lying on a surface ^ (cf. Figure B-1) . The tangent vector to the curve on the surface 
passing through P in the direction is given by 

s-.iL (B-u 

dX| 


/s 



Figure B-1 




ORIGINAL PAGS ^ 
OF POOR QUALITY 


The unit tangent vector is 


f. .iC 

ds, 


(B-2) 


vhere is the arc length in the direction. But from the chain rule, we obtain 


dr dr" dsi A 
TT- » 'T-' -rr^ • U 


dS| _ 

dXj ' dsj d X| " 'I dXj ' 

By definition = h^, the metric, so that 

dit 


t « 


h[ dXj 


(no sura on i) 


(B-3) 


The curve passing through point P in the direction e^ has a curvature 
associated with it which is given as 


K 


A A A 

dtj dtj dXj I dt| 

' ' dSj dXj dSj hj dXj 


I d_ / J_ _dT \ 

hi dfj \ hj dX[ / 




(B-4) 


d^f 


hj dXj dXj h| dXj 


The unit principal normal to the curve point at P is in the direction of and 
denoted by n^, while the unit normal to the surface N, at point P, is given by the 
cross product of the two tangent vectors at point P 


A A A 

t,x tg- Nsin«|/ 


(B-5) 


where tjj is the angle between t^^^ and t2» i.e. 


so that 


sin 


I ^ r I d r I dr 

cos y = t, ’ Ij, = -r X Z 1 — 

' 2 h, dX| hg dXg 

i/r = - cos^ ^ ^22 ^12 1 


(B-6) 
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where metric coefficients, 

and 


h,hj ■ 


•/^ 


II '•22 


The vector normal to both e^ and K is called the blnormal, B^, Hence the curvature 
can be represented by two components, one in the direction of N , and one in the 
direction of (cf. Figure B-2) . These components are called the normal curvature 


K^, and the geodesic curvature respectively 

= 1< • N * Kcos & 


Kg = K • B = K sin 6 


(7) 


where 0 is the angle between n and N. 



A A 

B = e, X N 


of page 


DESCRIPTION OF NACA FOUR-DIGIT AIRFOIL 

Consider an inertial reference frame in Cartesian coordinates 

= (xl, x^, x3) or equivalently (x, y, z) 

Attached to the airfoil is another Cartesian coordinate system given by 

= (x^, x^, x^) or (x, y, z) 

(cf. Figure B-3) . 
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The chord LT has length 

c = ((Xt - + (Xy 

and in general can be a function of z the spanwise coordinate, 

The thickness distribution for a MCA four-digit airfoil is given by 




+ Q| X + 


QgX'^ + OjX' 


+ Q 


a 

.X 


(B-8) 


(B-9) 


where x“x/c (cf kef. 50), For a cambered airfoil the thickness distribution is 
added onto the mean camber line. In the following, we will assume a symmetric 
uncambered airfoil with the mean line lying on the x axis. 

Hence the designation of the airfoil which is considered is NACA OOXX, 
where XX refers to the thickness ratio t/c. In order to obtain a single valued 
7 function, the polar angle 6 will be introduced (cf. Figure B-5). 



Figure B-5 
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The following relationships 

hold 





X = 

-C/2 

COS(^ d 

’C/2 

= -^(1- 
2 ^ 

• COS (p) 

rw 

X = 

JL 

C 

= i( 

1 - COS (ft) 


<}> ■ 0 leading edge 



O 
1 A 

(j) ~ r 

upper surface 

^ w ft trailing edge 



< 

7T - 

(j) - 2 tt 

lower surface 

Hence, Equation (9) becomes 







Jl 

c 

r*J 

r y = 

1 /c 

.20 




where g((}i) is the polar representation of the term in the brackets in Equation ( 

^ A 

and t and c can be functions of z. Defining g(^) ■ g(<}')/.2 and t(z) *• t(z)/c( 2 ) 
y^. becomes 

y = t{z)g((^) 

The coordinate x can also be represented as a function of z and ij), i.e. 

X = { ! - COS (jf>) = c(z)f(<^) 

and conversely 4' can be represented as a function of x and z 

The relationship between the (x,y) and (x,y) coordinate systems is 

"c" ” { (^~\)cosa ~ (y-yi_)sin a} 

— = ""^L^sina + (7-7^)005 0 } 

alternatively 

■q" " 

1^7^ + y COS a - X sin a I 

Given a radius vector from the origin of the (x,y,z) frame, 

r = X i + y j + zk 


(B-lO) 


(B-ll) 

-9) 

(B-12) 


(B-13) 

(B-14) 
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for the particular system under consideration 

x,+ c(5)f{(^)cosa + }(z)g{<j/>)sina 

i-J. 


(B-15) 


Hence 


■y + t(z)g((^)cosa - c(2) f (<^) slna 

1 =T(z,(j|!>) 

I ='z 


y =y(z,c/>) 

and the surface coordinates will be lines of constant and constant z. 


(B-16) 


In general on the surface we can write in symbolic form 

T + y z)J +T k (B-i7) 

The tangent vectors along the lines <(• “ constant and z >= constant can be computed 
by differentiating the radius vector with respect to the coordinate line, i. e. 



dT 

dx 0 ^ 



d(p 

d<p d(j3 


dr" 

(3x A , d y 4 

^z " 


2 Wl.l 1 I II ■ 1 

dz 

dz dz 

The unit tangent vectors are 

A 

T, 

X(^i + y«^] 

|t| 



1 z 

^ ft A 

Xj i + Vjj +k 

^z 

Xg i" y^^ 


The metric coefficients become 

= 8„ H-y/ 

~ ^ + Iz^ 


(B-18) 


(B-19) 


(B-20) 
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and 


9 12’ ^(jb^Z 


If ^2 orthogonal , S 0 

The normal to the surface Is obtained from the cross product of the surface unit 
tangent vectors 

A A A 

t| X “ N sin xj/ 

where <|/ is the angle between t^^ and 

In order to obtain sini;», the dot product of the unit tangent vectors is employed 

A A 

t| ■ tg = COS »// 


or 


l [xj? + yj4k]. i[,^t + y^f] = + 


h,l'2 


yielding the following relationship 
Employing the expression 


sin i// = - cos^^ 

and with some algebraic manipulation the desired result Is obtained. 

I r ->''2 

SI n r — — 'I g g - g ^ r 
h,hgV^M'^i2 ^12 J 

In the case where <{) and z are orthogonal to each other, ijj => % and g « 0 or 

z c>2^2 * 

= 0 

or 

Xj = = 0 

The normal to the surface can now be computed 

I r~A A _ __,at 

~ y<^ 


(B-21) 


A 

N = 


(B-22) 
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From Equations (D-4) and (B-20), the curvature can be obtained in the ^ direction 
(1 dlrecticn) 


J- iiiL .}JL + JL £L 

h|® d<p dcfi h|^ d(j3^ 



d 

dcf} 


[if/ + V] 


1/2 


d r An 

^[si+Ji+Jk] 


. ' r ? ^ 

“ ‘hj^ [v' -^7(^1] hV + y<^c/,n 


Similarly K 2 is given as 


T = _ „i_ dr* j_ 

^ hg^ dz dz h| dz^ 

~ [7^x22 yryzr][K 2 i -Hyzl + l!] + [%i+yzrj] (B-2/0 


The geodesic curvatures can now be computed 

Kg =!<, ■'b, = - K| ‘ (t, XN) 

where t^, and N are given by Equations (B-18) , (B-22) , and (B-23) . 

After some algebra, we obtain 

' *h, ~^7y(^) (B-25) 

Similarly K becomes 
§2 

" hi^^hg^'sin l// ~ "^2 y^JT^ ^21 "^yc^ yi2 (B-26) 


66 




We will now consider a NAGA 4 digit airfoil and further assure Xj^ is a constant 
(non swept trailing edge) . From Equation (15) we can obtain the appropriate 
derivatives 

X « + c(f) f (<^) cos a + t(z)g((/))sina 

“ c(x)f^cosa + tg^sina 

(B-27) 

J.~ - Cjfcosa + t^gsina 
x^z ' ^zz ^zz 


» ccosaf^^+ tsinag^^ 


y = + t(z)g{^)cosa - c(z)f(c^)sina 

^ = g^tc.osa - f^csina 

= g,^<^tcosa - f^<^csina (B-28) 

y- = g cosa tj " f sinac^ 
y__ = gcosatjj “ fsinaCg- 


For a rectangular plan form c and t are constant so that 



Substituting Equations (B-27) and (B-28) into Equations (B-18) to (B-24) the appropriate 
geometric quantities can be obtained specialized for a NACA OOXX airfoil. 
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For the present code, the curvature terms are not used directly, 
metric coefficients, the Jacobian, ,T, and the Christoffel sy^'^ls, 

are used. These geometric functions are presented below. 


Instead, the 
(cf Ref. 4) 


Metric Coefficients: 


Basis vectors 


g,i • h,* ■ s/ + 5/ 

955' ' 

'i • Y 

. x,i + +(< 

S ■ i {- Y + V “ 


W » x^Yj - y,^x- 


(B-29) 


(B-30) 


and 


e, = F|/h, , eg = eg/hg , 63 ‘ 63 


Christoffell Symbols (27 components) 


= i''^ = r^ = r^ = o 

^3 *23 31 32 

r ' = r ^ =0 

* 33 33 




1 = 1 , 2 , 3 


(B-31) 


(B-32a) 
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^12 = ^2i' ‘ i2{‘^^y2 ^ V<^z + y^y<^2^} 

Fg2 ■*■ y«^yzz^ ^z^4>^^^z\z ~ yzyfz^} 

^ » 

^13 ■' ^3] ' 73 Vj 

'’ 23 ' ^ 3 ! ■ 75'{[“f‘Vi * Vz’ ■“ '^^hH * ^hA + »+5^2>} 

^11 " ^^(f>^(p(p ” ^<p^<ptp^ 

^12 ’ ^21 “ j2 ^*«j!>y</)Z “ y^*(jbz ^ ^ 

-^(y-j -y^x._) 

r;l • ^3? ‘ 75 {“(y^y^^+ ‘ 

= '’3^ 73 {“(*^ Vz- y.#.y.(,z ) - “z <*/ ^^’>] 

r,? = -j ^ 

r,| = ^z, ■ -J^^^z-^^z'l 

^22 ~^V*zz ~ ^i>^ZZ^ 

These terms can be specialized for a NACA four digit airfoil by substituting 
Equations (B-27) and (B-28) into Equations (B-29) - (B-32) . 




-32 b) 
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Generalized Operator Compact Implicit Schemes 

In this section a procedure for generating generalized OCI schemes is reviewed 

Given 

L(u) ~ + b(x + c(x)u 

an expression relating L(u) and u is sought in the form 


-y- RUj = Q(Lu)j 




(C-l) 


where t. is the truncation error and Q and R are tridiagonal displacement 

J t 

operators. The maximum accuracy attainable is fourth order, i.e., t. O(h^). 

Expanding Eq. (20) in terms of and we obtain, 

J 3 . 


- [q‘ ( Lu)j_, + qJ(Lu)| + qJ'(Lu)j^,] 


(C-2) 


A Taylor series expansion yields for t . 


Tj - T° u(Xj) + Tj'u<”(Xj) + Tj^u‘^\xj) + Yu‘^’(Xj) 

+ Tj'‘u*''’(Xj) Xj) + T®u‘®’(xp + 0(h®) 


(C-3) 
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where superscripts in parentheses denote derivatives with respect to >:, and 
where 

I r 

Tf ' -p'll'-j' + fl' + fj> +h'(qfcj.|+ qjcj + qjcj,|)j 

T,‘ = 4“ [<'■)* ' '^l-l + =1) - h^qjcj^l -qJCj.,)] 


t/ ■ t''’’)" +>■]“)■< 'If + + <1/) 


(c-i) 


-h(q,*bj*i - qj'bj.i) - — (q|C)ti+ q]cj-|) 


-r*' 

Ti ' n 


-^(r; + (-D^p - +H)‘'''qfbj.|) 


■ (7^ ‘"f + (-l)"qj-) + b^(q;c,„ + (-l^qp,.,) 


1/ = 3, 4, 5, 6 


0 12 

For second order central finite differences we set T « T = T =0. This 


yields, when q? = 1 and q. = q. =0, the following relations 
1 11 


c 

^1 ■ 


-(rj‘ + rp + h^(q[Cj.| + qpj + qpj*,) 


rj^ - Tj" = hbj = Rcj 

r; + r,- = 2 


which recovers Eq. (45), i.e., 


r j“ = I - RCj / 2 
rj*' = I + Rcj/2 
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To obtain the fourth order operator compact implicit scheme we again set 

0 12 - c + 

X «=T «T *>0to obtain three expressions for r.’ ' in terms of the 

•• 0 

, (note that q , q f 0 and q is not necessarily unity), i.e., 

rj “ + rp + h2(qj"C|.| + qj^Cj +q+cj^,) (c-5a) 

rl - r- = Rcj_,qj" + Rc^qP + Rcj+iq^"^ + h^(qpj^, + qj" Cj.,) (c-5b) 

rj^ + rj" = 2(qj“+q5 +qp +2[Rcj^,q+ - Rcj.,qj-] + h2(qpj+, 4- qj"Cj.,) (C-5c) 

3 A c ~f* 

Now T and T must be set to evaluate q.’ ’ . The standard Swartz OCI 

-1 4 J 

method requires T = T =0, 

-i (r*-r-)- -l[Rcj^,q; + Rcj.,q,-] 

-h2(q/cj*|-qj-Cj.|) 


- [q/-qj-j =0 

(C-5d) 


(C-5e) 


Substi- 


-^{r* +r-)--i-[R>q/ - Rcj_,qf] - 4" [ 'if + <fr] = ° 

- h2(qp.^, + qfcj_,) 

rsj ( 5 ) C 6 ) \ A 

and results in a leading truncation error of the form (Au + Bu )h , 

■f - 

tuting (C-5b) and (C-5c) into (C--5d) and (C-53) r. and r. can be eliminated and a 

— c + 

system of two equations in q. and q. with q. as a parameter is obtained. The 

c J d 3 

parameter q^ is proportional to the determinant of the system. The values of 

q.’*^’"^ and are presented in Table I. 

As shown in Ref. 20 a cell Reynolds number stability condition exists 
for the Swartz OCI scheme, i.e., for Rc - ./l2 nonrealistic or oscillatory 


solutions will be obtained. In order to eliminate this restriction one can 
Dnc 

4 


3 4 (3) (4) 

relax the conditions T = T = 0, and allow the coefficients of u and u 
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"* C *1* 

By expanding 'I j ’ ’ ^ series in Rc 

q-'"'* • I X1''-*RC,"" (c-6) 

m=o * 

^ c 

12 parameters, Ajj,’ ’ m « 0, 1, 2, 3, are introduced. The equations for 

“ O(h^) and = O(h^) yield 5 linear relations, leaving at one's disposal 
"6 free” parameters plus a factor. 

These parameters can be set according to some criteria that would yield 
certain desirable properties for the difference equations. The following 
constraints are prescribed 


qj^ > 0, qj“ > 0, qj^ > 0 

q[bj +q/bj„ 

rj*> fj" >0 

~ r + r." 

1 J J 


(c-7) 


and h is sufficiently small so that 

lObj -bj.|-bj^i > 0 and 2 + hCj^, /bj^i > 0 
for J=» 2 , • •• ,J and Cj <o 

These conditions assure that R is diagonally dominant and Q is invertible for 
all R^. Further details are given in Ref. 27. The significance of this 
approach is that one can construct a scheme possessing certain desired proper- 
ties by employing a set of preassigned rules. This is contrary to usual practice, 
in which a scheme is chosen, and then its properties are determined. Although 
the computational effort in computing the q and r coefficients is not cheap, 
the actual percentage of the total cost is minimal. This has been borne out 
in actual computations. 

The q and r coefficients for the generalized OCI scheme described in 
Ref. 27 are given in Table II. 
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TABLE I, - OPERATOR COEFFICIEMTS FOR STANDARD 
OPERATOR COMPACT IMPLICIT SCHEME 

qj- ■ 6 - 5/>j + -PjPjt, 

qj' • 60 + I6/.J*, - 

q/ ■ 6 + - Zp^., - p^p^., 

■ "(fj* + '■f) + h®(qfCj.| + q'cj + qj*Cj^|) 

■ ‘'/■('"i'/’j-i) "TPi ) + '>)'■( 

where 

P\ “ 
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TABLE II. - OPERATOR COEFFZCIEKTS FOR GENERALIZED 
OPERATOR COmCT IMPLICIT SCHEME 

qf - 6 + [p,- 3 ]rcj +[pj]rcj* 
qj‘ • 60 +[|0 p,]ri:j + [pJrc' + [tj„p^]rc® 
qj‘ • 6 + [ p, + 3 ]rcj + [ p, + p J Rcj' + [ p J rcj’ 

where 

p, . 3 , Pj. - 0 , Pj ■ mox[7r, .iTj] 


’"i ■ + + '^J.|P| + ^1 ■ '5 -2Pj + (°'j-*)P| -3(T-j„ + crj) + 7T, 


‘W. 


0 

CTj >0 

_ 1 
TT_ ■ 

0 

2p^'~a^>0 


o^j <0 

Z 


2p,-cr2<0 


I n 


P^ " 


P./3 


[' + V] 


20*2 » 3tj_, - + 10 + 2hTj 


ffH-N 

-'Uj., I 

, ■ Pj-'^i + ’^l +2Tj.|(2 + h-^)pj 


with h sufficiently small so that 

IObj-bj.,-b^^, >0 and 2 + bcj^,/bj^, > O for J- 2 , • •• , J ond Cj^O 

where 

Rcj - hbj 

r7, Tj, r^ given in TABLE I 


rr = ROje-'"l/(l-e-'"‘) 
<■* = Rcj/{l-e-"'l) 

= - Rcj + cj 
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TABLE IV. - OI'ERATOR COEFFICIENTS FOR EL-MlSTIKy»VA’ VJERLE 
EXPONENTIAL BOX SCHEME 


rj“ - exp( -/J- >] 

rj'*^ « / [l‘-exp{ “^'*’)] 

Qj" = (I - rp /C2.p~) 

q; • /(2p^) 

qj “• qf + qf 
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TABLE V. - OPERATOR COEFFICIENTS FOR SECOND ORDER 

FINITE DIFFERENCES WITH ARTIFICIAL VISCOSITY 



whare S = Rc /2 for Rc, > Rc 

max j ' max 


^ " ~^^max'^^ Rc. < Rc^qj^ 


Rc 


max 


IRcJ 


and RCj = hbj 


foz' IRCjl < RCrna,, rj, r J , r| , qT, qj, qj 
reduce to standard finite differences. 

where Rc is the maximum allowable cell Reynolds number 
max ■' 
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IG. 3. NEWMAN AIRFOIL SERIES 2. SKIN FRICTION DISTRIBUTION 


AS A FUNCTION OF x 
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□ FOURTH ORDER G/OCI 
A SECOND ORDER F/ D 







LIGHTHILL LOW FREQUENCY 


LIGHTHILL HIGH FREQUENCY 



iO = UJX / u. 


FIG. 7. PHASE ANGLE BETWEEN EDGE VELOCITY AND WALL SHEAR STRESS 

AS A FUNCTION OF REDUCED FREQUENCY 
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